corentinchoisy/simul_these
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Generated N=50 simulated responses

Browse Source
This commit is contained in:
Corentin Choisy
2024-02-22 16:30:03 +01:00
parent 6e0f980484
commit da3fa74ff5
8 changed files with 13149 additions and 602 deletions
Download Patch File Download Diff File Expand all files Collapse all files

BIN
RProject/.RData
View File

Binary file not shown.

1010
RProject/.Rhistory
View File

File diff suppressed because it is too large Load Diff

196
RProject/Scripts/pcm.R
View File

@ -2266,6 +2266,8 @@ for (x in results[seq(2,length(results))]) {
res.dat.dif <- bind_rows(res.dat.dif,y)
}
res.dat$bias <- res.dat$eff.size-res.dat$m.beta
res.dat.dif$bias <- res.dat.dif$eff.size-res.dat.dif$m.beta
##############################################################################
#----------------------------------------------------------------------------#
@ -2281,118 +2283,118 @@ res.dat$theoretical.power <- 0
## Scénarios J=4 / M=2
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'A') & res.dat$N==100,]$theoretical.power <- 0.0250
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'B') & res.dat$N==100,]$theoretical.power <- 0.1524
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'C') & res.dat$N==100,]$theoretical.power <- 0.1524
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'D') & res.dat$N==100,]$theoretical.power <- 0.4626
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'E') & res.dat$N==100,]$theoretical.power <- 0.4626
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'F') & res.dat$N==100,]$theoretical.power <- 0.1524
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'G') & res.dat$N==100,]$theoretical.power <- 0.4626
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'A') & res.dat$N==100,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'B') & res.dat$N==100,]$theoretical.power <- 0.1543
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'C') & res.dat$N==100,]$theoretical.power <- 0.1543
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'D') & res.dat$N==100,]$theoretical.power <- 0.4627
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'E') & res.dat$N==100,]$theoretical.power <- 0.4627
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'F') & res.dat$N==100,]$theoretical.power <- 0.1543
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'G') & res.dat$N==100,]$theoretical.power <- 0.4627
res.dat[res.dat$scenario %in% paste0(1,'C') & res.dat$N==100,]$theoretical.power <- 0.4626
res.dat[res.dat$scenario %in% paste0(1,'D') & res.dat$N==100,]$theoretical.power <- 0.1524
res.dat[res.dat$scenario %in% paste0(1,'E') & res.dat$N==100,]$theoretical.power <- 0.4626
res.dat[res.dat$scenario %in% paste0(1,'C') & res.dat$N==100,]$theoretical.power <- 0.4627
res.dat[res.dat$scenario %in% paste0(1,'D') & res.dat$N==100,]$theoretical.power <- 0.1543
res.dat[res.dat$scenario %in% paste0(1,'E') & res.dat$N==100,]$theoretical.power <- 0.4627
## Scénarios J=4 / M=4
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'A') & res.dat$N==100,]$theoretical.power <- 0.0250
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'B') & res.dat$N==100,]$theoretical.power <- 0.1522
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'C') & res.dat$N==100,]$theoretical.power <- 0.1522
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'D') & res.dat$N==100,]$theoretical.power <- 0.4626
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'E') & res.dat$N==100,]$theoretical.power <- 0.4626
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'F') & res.dat$N==100,]$theoretical.power <- 0.1522
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'G') & res.dat$N==100,]$theoretical.power <- 0.4626
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'A') & res.dat$N==100,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'B') & res.dat$N==100,]$theoretical.power <- 0.2177
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'C') & res.dat$N==100,]$theoretical.power <- 0.2177
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'D') & res.dat$N==100,]$theoretical.power <- 0.6586
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'E') & res.dat$N==100,]$theoretical.power <- 0.6586
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'F') & res.dat$N==100,]$theoretical.power <- 0.2177
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'G') & res.dat$N==100,]$theoretical.power <- 0.6586
res.dat[res.dat$scenario %in% paste0(2,'C') & res.dat$N==100,]$theoretical.power <- 0.4626
res.dat[res.dat$scenario %in% paste0(2,'D') & res.dat$N==100,]$theoretical.power <- 0.1522
res.dat[res.dat$scenario %in% paste0(2,'E') & res.dat$N==100,]$theoretical.power <- 0.4626
res.dat[res.dat$scenario %in% paste0(2,'C') & res.dat$N==100,]$theoretical.power <- 0.6586
res.dat[res.dat$scenario %in% paste0(2,'D') & res.dat$N==100,]$theoretical.power <- 0.2177
res.dat[res.dat$scenario %in% paste0(2,'E') & res.dat$N==100,]$theoretical.power <- 0.6586
## Scénarios J=7 / M=2
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'A') & res.dat$N==100,]$theoretical.power <- 0.0250
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'B') & res.dat$N==100,]$theoretical.power <- 0.1857
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'C') & res.dat$N==100,]$theoretical.power <- 0.1857
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'D') & res.dat$N==100,]$theoretical.power <- 0.5669
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'E') & res.dat$N==100,]$theoretical.power <- 0.5669
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'F') & res.dat$N==100,]$theoretical.power <- 0.1857
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'G') & res.dat$N==100,]$theoretical.power <- 0.5669
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'A') & res.dat$N==100,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'B') & res.dat$N==100,]$theoretical.power <- 0.1870
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'C') & res.dat$N==100,]$theoretical.power <- 0.1870
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'D') & res.dat$N==100,]$theoretical.power <- 0.5666
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'E') & res.dat$N==100,]$theoretical.power <- 0.5666
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'F') & res.dat$N==100,]$theoretical.power <- 0.1870
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'G') & res.dat$N==100,]$theoretical.power <- 0.5666
res.dat[res.dat$scenario %in% paste0(3,'C') & res.dat$N==100,]$theoretical.power <- 0.5669
res.dat[res.dat$scenario %in% paste0(3,'D') & res.dat$N==100,]$theoretical.power <- 0.1857
res.dat[res.dat$scenario %in% paste0(3,'E') & res.dat$N==100,]$theoretical.power <- 0.5669
res.dat[res.dat$scenario %in% paste0(3,'C') & res.dat$N==100,]$theoretical.power <- 0.5666
res.dat[res.dat$scenario %in% paste0(3,'D') & res.dat$N==100,]$theoretical.power <- 0.1870
res.dat[res.dat$scenario %in% paste0(3,'E') & res.dat$N==100,]$theoretical.power <- 0.5666
## Scénarios J=7 / M=4
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'A') & res.dat$N==100,]$theoretical.power <- 0.0250
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'B') & res.dat$N==100,]$theoretical.power <- 0.1852
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'C') & res.dat$N==100,]$theoretical.power <- 0.1852
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'D') & res.dat$N==100,]$theoretical.power <- 0.5669
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'E') & res.dat$N==100,]$theoretical.power <- 0.5669
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'F') & res.dat$N==100,]$theoretical.power <- 0.1852
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'G') & res.dat$N==100,]$theoretical.power <- 0.5669
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'A') & res.dat$N==100,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'B') & res.dat$N==100,]$theoretical.power <- 0.2450
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'C') & res.dat$N==100,]$theoretical.power <- 0.2450
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'D') & res.dat$N==100,]$theoretical.power <- 0.7136
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'E') & res.dat$N==100,]$theoretical.power <- 0.7136
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'F') & res.dat$N==100,]$theoretical.power <- 0.2450
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'G') & res.dat$N==100,]$theoretical.power <- 0.7136
res.dat[res.dat$scenario %in% paste0(4,'C') & res.dat$N==100,]$theoretical.power <- 0.5669
res.dat[res.dat$scenario %in% paste0(4,'D') & res.dat$N==100,]$theoretical.power <- 0.1852
res.dat[res.dat$scenario %in% paste0(4,'E') & res.dat$N==100,]$theoretical.power <- 0.5669
res.dat[res.dat$scenario %in% paste0(4,'C') & res.dat$N==100,]$theoretical.power <- 0.7136
res.dat[res.dat$scenario %in% paste0(4,'D') & res.dat$N==100,]$theoretical.power <- 0.2450
res.dat[res.dat$scenario %in% paste0(4,'E') & res.dat$N==100,]$theoretical.power <- 0.7136
### Scénarios N=200
## Scénarios J=4 / M=2
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'A') & res.dat$N==200,]$theoretical.power <- 0.0250
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'B') & res.dat$N==200,]$theoretical.power <- 0.2613
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'C') & res.dat$N==200,]$theoretical.power <- 0.2613
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'A') & res.dat$N==200,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'B') & res.dat$N==200,]$theoretical.power <- 0.2618
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'C') & res.dat$N==200,]$theoretical.power <- 0.2618
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'D') & res.dat$N==200,]$theoretical.power <- 0.7507
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'E') & res.dat$N==200,]$theoretical.power <- 0.7507
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'F') & res.dat$N==200,]$theoretical.power <- 0.2613
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'F') & res.dat$N==200,]$theoretical.power <- 0.2618
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'G') & res.dat$N==200,]$theoretical.power <- 0.7507
res.dat[res.dat$scenario %in% paste0(1,'C') & res.dat$N==200,]$theoretical.power <- 0.7507
res.dat[res.dat$scenario %in% paste0(1,'D') & res.dat$N==200,]$theoretical.power <- 0.2613
res.dat[res.dat$scenario %in% paste0(1,'D') & res.dat$N==200,]$theoretical.power <- 0.2618
res.dat[res.dat$scenario %in% paste0(1,'E') & res.dat$N==200,]$theoretical.power <- 0.7507
## Scénarios J=4 / M=4
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'A') & res.dat$N==200,]$theoretical.power <- 0.0250
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'B') & res.dat$N==200,]$theoretical.power <- 0.2606
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'C') & res.dat$N==200,]$theoretical.power <- 0.2606
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'D') & res.dat$N==200,]$theoretical.power <- 0.7507
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'E') & res.dat$N==200,]$theoretical.power <- 0.7507
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'F') & res.dat$N==200,]$theoretical.power <- 0.2606
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'G') & res.dat$N==200,]$theoretical.power <- 0.4626
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'A') & res.dat$N==200,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'B') & res.dat$N==200,]$theoretical.power <- 0.3875
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'C') & res.dat$N==200,]$theoretical.power <- 0.3875
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'D') & res.dat$N==200,]$theoretical.power <- 0.9161
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'E') & res.dat$N==200,]$theoretical.power <- 0.9161
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'F') & res.dat$N==200,]$theoretical.power <- 0.3875
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'G') & res.dat$N==200,]$theoretical.power <- 0.9161
res.dat[res.dat$scenario %in% paste0(2,'C') & res.dat$N==200,]$theoretical.power <- 0.7507
res.dat[res.dat$scenario %in% paste0(2,'D') & res.dat$N==200,]$theoretical.power <- 0.2606
res.dat[res.dat$scenario %in% paste0(2,'E') & res.dat$N==200,]$theoretical.power <- 0.7507
res.dat[res.dat$scenario %in% paste0(2,'C') & res.dat$N==200,]$theoretical.power <- 0.9161
res.dat[res.dat$scenario %in% paste0(2,'D') & res.dat$N==200,]$theoretical.power <- 0.3875
res.dat[res.dat$scenario %in% paste0(2,'E') & res.dat$N==200,]$theoretical.power <- 0.9161
## Scénarios J=7 / M=2
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'A') & res.dat$N==200,]$theoretical.power <- 0.0250
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'B') & res.dat$N==200,]$theoretical.power <- 0.3256
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'C') & res.dat$N==200,]$theoretical.power <- 0.3256
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'A') & res.dat$N==200,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'B') & res.dat$N==200,]$theoretical.power <- 0.3258
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'C') & res.dat$N==200,]$theoretical.power <- 0.3258
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'D') & res.dat$N==200,]$theoretical.power <- 0.8538
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'E') & res.dat$N==200,]$theoretical.power <- 0.8538
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'F') & res.dat$N==200,]$theoretical.power <- 0.3256
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'F') & res.dat$N==200,]$theoretical.power <- 0.3258
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'G') & res.dat$N==200,]$theoretical.power <- 0.8538
res.dat[res.dat$scenario %in% paste0(3,'C') & res.dat$N==200,]$theoretical.power <- 0.8538
res.dat[res.dat$scenario %in% paste0(3,'D') & res.dat$N==200,]$theoretical.power <- 0.3256
res.dat[res.dat$scenario %in% paste0(3,'D') & res.dat$N==200,]$theoretical.power <- 0.3258
res.dat[res.dat$scenario %in% paste0(3,'E') & res.dat$N==200,]$theoretical.power <- 0.8538
## Scénarios J=7 / M=4
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'A') & res.dat$N==200,]$theoretical.power <- 0.0250
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'B') & res.dat$N==200,]$theoretical.power <- 0.3251
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'C') & res.dat$N==200,]$theoretical.power <- 0.3251
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'D') & res.dat$N==200,]$theoretical.power <- 0.8538
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'E') & res.dat$N==200,]$theoretical.power <- 0.8538
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'F') & res.dat$N==200,]$theoretical.power <- 0.3251
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'G') & res.dat$N==200,]$theoretical.power <- 0.8538
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'A') & res.dat$N==200,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'B') & res.dat$N==200,]$theoretical.power <- 0.4321
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'C') & res.dat$N==200,]$theoretical.power <- 0.4321
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'D') & res.dat$N==200,]$theoretical.power <- 0.9471
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'E') & res.dat$N==200,]$theoretical.power <- 0.9471
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'F') & res.dat$N==200,]$theoretical.power <- 0.4321
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'G') & res.dat$N==200,]$theoretical.power <- 0.9471
res.dat[res.dat$scenario %in% paste0(4,'C') & res.dat$N==200,]$theoretical.power <- 0.8538
res.dat[res.dat$scenario %in% paste0(4,'D') & res.dat$N==200,]$theoretical.power <- 0.3251
res.dat[res.dat$scenario %in% paste0(4,'E') & res.dat$N==200,]$theoretical.power <- 0.8538
res.dat[res.dat$scenario %in% paste0(4,'C') & res.dat$N==200,]$theoretical.power <- 0.9471
res.dat[res.dat$scenario %in% paste0(4,'D') & res.dat$N==200,]$theoretical.power <- 0.4321
res.dat[res.dat$scenario %in% paste0(4,'E') & res.dat$N==200,]$theoretical.power <- 0.9471
@ -2401,35 +2403,35 @@ res.dat[res.dat$scenario %in% paste0(4,'E') & res.dat$N==200,]$theoretical.power
## Scénarios J=4 / M=2
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'A') & res.dat$N==300,]$theoretical.power <- 0.0250
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'B') & res.dat$N==300,]$theoretical.power <- 0.3658
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'C') & res.dat$N==300,]$theoretical.power <- 0.3658
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'A') & res.dat$N==300,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(1,5,7,9,11),'B') & res.dat$N==300,]$theoretical.power <- 0.3660
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'C') & res.dat$N==300,]$theoretical.power <- 0.3660
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'D') & res.dat$N==300,]$theoretical.power <- 0.8981
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'E') & res.dat$N==300,]$theoretical.power <- 0.8981
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'F') & res.dat$N==300,]$theoretical.power <- 0.3658
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'F') & res.dat$N==300,]$theoretical.power <- 0.3660
res.dat[res.dat$scenario %in% paste0(c(5,7,9,11),'G') & res.dat$N==300,]$theoretical.power <- 0.8981
res.dat[res.dat$scenario %in% paste0(1,'C') & res.dat$N==300,]$theoretical.power <- 0.8981
res.dat[res.dat$scenario %in% paste0(1,'D') & res.dat$N==300,]$theoretical.power <- 0.3658
res.dat[res.dat$scenario %in% paste0(1,'D') & res.dat$N==300,]$theoretical.power <- 0.3660
res.dat[res.dat$scenario %in% paste0(1,'E') & res.dat$N==300,]$theoretical.power <- 0.8981
## Scénarios J=4 / M=4
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'A') & res.dat$N==300,]$theoretical.power <- 0.0250
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'B') & res.dat$N==300,]$theoretical.power <- 0.3652
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'C') & res.dat$N==300,]$theoretical.power <- 0.3652
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'D') & res.dat$N==300,]$theoretical.power <- 0.8981
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'E') & res.dat$N==300,]$theoretical.power <- 0.8981
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'F') & res.dat$N==300,]$theoretical.power <- 0.3652
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'G') & res.dat$N==300,]$theoretical.power <- 0.8981
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'A') & res.dat$N==300,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(2,6,8,10,12),'B') & res.dat$N==300,]$theoretical.power <- 0.5373
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'C') & res.dat$N==300,]$theoretical.power <- 0.5373
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'D') & res.dat$N==300,]$theoretical.power <- 0.9834
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'E') & res.dat$N==300,]$theoretical.power <- 0.9834
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'F') & res.dat$N==300,]$theoretical.power <- 0.5373
res.dat[res.dat$scenario %in% paste0(c(6,8,10,12),'G') & res.dat$N==300,]$theoretical.power <- 0.9834
res.dat[res.dat$scenario %in% paste0(2,'C') & res.dat$N==300,]$theoretical.power <- 0.8981
res.dat[res.dat$scenario %in% paste0(2,'D') & res.dat$N==300,]$theoretical.power <- 0.3652
res.dat[res.dat$scenario %in% paste0(2,'E') & res.dat$N==300,]$theoretical.power <- 0.8981
res.dat[res.dat$scenario %in% paste0(2,'C') & res.dat$N==300,]$theoretical.power <- 0.9834
res.dat[res.dat$scenario %in% paste0(2,'D') & res.dat$N==300,]$theoretical.power <- 0.5373
res.dat[res.dat$scenario %in% paste0(2,'E') & res.dat$N==300,]$theoretical.power <- 0.9834
## Scénarios J=7 / M=2
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'A') & res.dat$N==300,]$theoretical.power <- 0.0250
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'A') & res.dat$N==300,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(3,13,15,17,19),'B') & res.dat$N==300,]$theoretical.power <- 0.4550
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'C') & res.dat$N==300,]$theoretical.power <- 0.4550
res.dat[res.dat$scenario %in% paste0(c(13,15,17,19),'D') & res.dat$N==300,]$theoretical.power <- 0.9584
@ -2443,19 +2445,19 @@ res.dat[res.dat$scenario %in% paste0(3,'E') & res.dat$N==300,]$theoretical.power
## Scénarios J=7 / M=4
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'A') & res.dat$N==300,]$theoretical.power <- 0.0250
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'B') & res.dat$N==300,]$theoretical.power <- 0.4547
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'C') & res.dat$N==300,]$theoretical.power <- 0.4547
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'D') & res.dat$N==300,]$theoretical.power <- 0.9584
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'E') & res.dat$N==300,]$theoretical.power <- 0.9584
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'F') & res.dat$N==300,]$theoretical.power <- 0.4547
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'G') & res.dat$N==300,]$theoretical.power <- 0.9584
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'A') & res.dat$N==300,]$theoretical.power <- 0.05
res.dat[res.dat$scenario %in% paste0(c(4,14,16,18,20),'B') & res.dat$N==300,]$theoretical.power <- 0.5907
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'C') & res.dat$N==300,]$theoretical.power <- 0.5907
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'D') & res.dat$N==300,]$theoretical.power <- 0.9919
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'E') & res.dat$N==300,]$theoretical.power <- 0.9919
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'F') & res.dat$N==300,]$theoretical.power <- 0.5907
res.dat[res.dat$scenario %in% paste0(c(14,16,18,20),'G') & res.dat$N==300,]$theoretical.power <- 0.9919
res.dat[res.dat$scenario %in% paste0(4,'C') & res.dat$N==300,]$theoretical.power <- 0.9584
res.dat[res.dat$scenario %in% paste0(4,'D') & res.dat$N==300,]$theoretical.power <- 0.4547
res.dat[res.dat$scenario %in% paste0(4,'E') & res.dat$N==300,]$theoretical.power <- 0.9584
res.dat[res.dat$scenario %in% paste0(4,'C') & res.dat$N==300,]$theoretical.power <- 0.9919
res.dat[res.dat$scenario %in% paste0(4,'D') & res.dat$N==300,]$theoretical.power <- 0.5907
res.dat[res.dat$scenario %in% paste0(4,'E') & res.dat$N==300,]$theoretical.power <- 0.9919
### DIF scenarios
res.dat.dif$theoretical.power <- res.dat[61:nrow(res.dat),]$theoretical.power
res.dat.dif$theoretical.power <- res.dat[61:nrow(res.dat),]$theoretical.power

293
RProject/Scripts/power_analysis.R Normal file
View File

@ -0,0 +1,293 @@
##############################################################################
#----------------------------------------------------------------------------#
################################# POWER PLOTS ################################
#----------------------------------------------------------------------------#
##############################################################################
## Plot 1 - baseline scenarios vs theoretical power
par(mfrow=c(1,2))
# theoretical
plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),rep(unique(plot.dat[plot.dat$scenario.type=="A",]$theoretical.power),3),
xlab='N',ylab='Theoretical power',
xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(0,1),lty=4,col='#03a18a')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c6d18d',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#da77c7',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#b5a180',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#30a466',lty=4)
lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',lty=4)
points(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$N),
rep(unique(plot.dat[plot.dat$scenario %in% c("1A",'1A'),]$theoretical.power),3),col='#03a18a',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$theoretical.power),col='#c6d18d',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$theoretical.power),col='#c0c23b',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$theoretical.power),col='#da77c7',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$theoretical.power),col='#a12471',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$theoretical.power),col='#b5a180',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$theoretical.power),col='#9b6541',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$theoretical.power),col='#30a466',pch=17)
points(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$theoretical.power),col='#1a342b',pch=17)
legend("topleft",lty=rep(4,5),col=c("#03a18a","#c0c23b","#a12471","#9b6541","#1a342b"),
pch=rep(17,5),
legend=c("Scenario A",
"Scenario 1-2 / B-D",
"Scenario 3-4 / B-D",
"Scenario 1-2 / C-E",
"Scenario 3-4 / C-E"),cex=0.7)
# real
plot.dat <- res.dat[res.dat$scenario %in% sapply(c('A','B','C','D','E'),function(x) paste0(1:4,x)),]
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
xlab='N',ylab='Proportion of null hypothesis rejection',
xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(0,1),lty=4,col='#03a18a')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
mean.A,col='#c6d18d',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
mean.A,col='#c0c23b',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
mean.A,col='#da77c7',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
mean.A,col='#a12471',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
mean.A,col='#b5a180',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
mean.A,col='#9b6541',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
mean.A,col='#30a466',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
mean.A,col='#1a342b',lty=4)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario.type=="A" & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,pch=17,col='#03a18a')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1B","1D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("1B",'1D'),]$N),
mean.A,col='#c6d18d',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2B","2D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("2B",'2D'),]$N),
mean.A,col='#c0c23b',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3B","3D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("3B",'3D'),]$N),
mean.A,col='#da77c7',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4B","4D") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("4B",'4D'),]$N),
mean.A,col='#a12471',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("1C","1E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("1C",'1E'),]$N),
mean.A,col='#b5a180',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("2C","2E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("2C",'2E'),]$N),
mean.A,col='#9b6541',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("3C","3E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("3C",'3E'),]$N),
mean.A,col='#30a466',pch=17)
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario%in% c("4C","4E") & plot.dat$N==x,]$h0.rejected.p))
points(unique(plot.dat[plot.dat$scenario %in% c("4C",'4E'),]$N),
mean.A,col='#1a342b',pch=17)
legend("topleft",lty=rep(4,5),col=c("#03a18a","#c0c23b","#a12471","#9b6541","#1a342b"),
pch=rep(17,5),
legend=c("Scenario A",
"Scenario 1-2 / B-D",
"Scenario 3-4 / B-D",
"Scenario 1-2 / C-E",
"Scenario 3-4 / C-E"),cex=0.7)
##############################################################################
#----------------------------------------------------------------------------#
################################# ALPHA PLOTS ################################
#----------------------------------------------------------------------------#
##############################################################################
par(mfrow=c(1,1))
#### Scenarios with J=4 / M=2
# A / H=0
plot.dat <- res.dat
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="1A" & plot.dat$N==x,]$h0.rejected.p))
plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
xlab='N',ylab='Proportion of null hypothesis rejection',
xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(0,1),lty=4,col='#03a18a',main='Scenarios where H0 is TRUE')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#03a18a')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="5A" & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#a12471')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#a12471')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="7A" & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#9b6541')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#9b6541')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="9A" & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#1a342b')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#1a342b')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="11A" & plot.dat$N==x,]$h0.rejected.p))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#c0c23b')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#c0c23b')
##############################################################################
#----------------------------------------------------------------------------#
################################# BIAS PLOTS #################################
#----------------------------------------------------------------------------#
##############################################################################
par(mfrow=c(1,1))
#### Scenarios with J=4 / M=2
# A / H=0
plot.dat <- res.dat
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="1A" & plot.dat$N==x,]$bias))
plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
xlab='N',ylab='Bias',
xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(-1,1),lty=4,col='#03a18a',main='Scenarios where H0 is TRUE')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#03a18a')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="5A" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#a12471')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#a12471')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="7A" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#9b6541')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#9b6541')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="9A" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#1a342b')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#1a342b')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="11A" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#c0c23b')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#c0c23b')
par(mfrow=c(1,1))
#### Scenarios with J=4 / M=2
# B/ H=0
plot.dat <- res.dat
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="1B" & plot.dat$N==x,]$bias))
plot(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
xlab='N',ylab='Bias',
xaxt='n',yaxt='n',type='l',xlim=c(100,300),ylim=c(-1,1),lty=4,col='#03a18a',main='Scenarios where H0 is TRUE')
axis(1,c(100,200,300))
axis(2,seq(0,1,0.1))
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#03a18a')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="5B" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#a12471')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#a12471')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="7B" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#9b6541')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#9b6541')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="9B" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#1a342b')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#1a342b')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="11B" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='#c0c23b')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='#c0c23b')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="5C" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='red')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='red')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="7C" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='blue')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='blue')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="9C" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='purple')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='purple')
mean.A <- sapply(c(100,200,300),function(x) mean(plot.dat[plot.dat$scenario=="11C" & plot.dat$N==x,]$bias))
lines(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,lty=4,col='orange')
points(unique(plot.dat[plot.dat$scenario.type=="A",]$N),mean.A,
pch=17,col='orange')

1515
Scripts/Analysis/DIF/pcm_dif_50.do Normal file
View File

File diff suppressed because it is too large Load Diff

3815
Scripts/Scenarios/DIF/scenarios_DIF_J4_baseline_50.do Normal file
View File

File diff suppressed because it is too large Load Diff

6109
Scripts/Scenarios/DIF/scenarios_DIF_J7_baseline_50.do Normal file
View File

File diff suppressed because it is too large Load Diff

813
Scripts/Scenarios/NoDIF/scenarios_noDIF_baseline_50.do Normal file
View File

@ -0,0 +1,813 @@
*=================================================================================================================================================
* Date : 2024-01-04
* Stata version : Stata 18 SE
*
* This program creates dataset without DIF for a randomized controlled trial scenario
*
* ado-files needed : - simirt (version 4.3 August 29, 2019, available on OSF)
*
* outputs : scenario_1,scenario_2,scenario_3,scenario_4, for N=50/200/300
*
*
* Warning : To obtain reproduce the data obtained in the .csv files in this repository, use 'simirt_setseed.ado' instead of 'simirt.ado'
*
*
*================================================================================================================================================
* Load simirt.ado
adopath+"/home/corentin/Documents/These/Recherche/Simulations/Modules/"
* Set data output folder path
local path = "/home/corentin/Documents/These/Recherche/Simulations/Data/NoDIF"
* Scenarios with : n = 100
*==========================
local path = "/home/corentin/Documents/These/Recherche/Simulations/Data/NoDIF/N50"
local Nn = 50
** Scenario 1: J = 4 items / M = 2 modalities / N=50 per group / TT=treatment variable
* Scenario 1A : H_0 is TRUE
di "SCENARIO 1A - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-0.84 \ -0.25 \ 0.25 \ 0.84)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(4) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-0.84 \ -0.25 \ 0.25 \ 0.84)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(4) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_1A_50.csv", replace
* Scenario 1B : H_0 is FALSE / Effect size = 0.2
di "SCENARIO 1B - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-0.84 \ -0.25 \ 0.25 \ 0.84)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(4) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-0.84 \ -0.25 \ 0.25 \ 0.84)
qui simirt, nbobs(`Nn') mu(0.2) cov(1) dim(4) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_1B_50.csv", replace
* Scenario 1C : H_0 is FALSE / Effect size = 0.4
di "SCENARIO 1C - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-0.84 \ -0.25 \ 0.25 \ 0.84)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(4) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-0.84 \ -0.25 \ 0.25 \ 0.84)
qui simirt, nbobs(`Nn') mu(0.4) cov(1) dim(4) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_1C_50.csv", replace
* Scenario 1D : H_0 is FALSE / Effect size = -0.2
di "SCENARIO 1D - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-0.84 \ -0.25 \ 0.25 \ 0.84)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(4) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-0.84 \ -0.25 \ 0.25 \ 0.84)
qui simirt, nbobs(`Nn') mu(-0.2) cov(1) dim(4) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_1D_50.csv", replace
* Scenario 1E : H_0 is FALSE / Effect size = -0.4
di "SCENARIO 1E - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-0.84 \ -0.25 \ 0.25 \ 0.84)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(4) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-0.84 \ -0.25 \ 0.25 \ 0.84)
qui simirt, nbobs(`Nn') mu(-0.4) cov(1) dim(4) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_1E_50.csv", replace
**------------------------------------------------------------------------------------**
** Scenario 2: J = 4 items / M = 4 modalities / N=50 per group / TT=treatment variable
* Scenario 2A : H_0 is TRUE
di "SCENARIO 2A - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-1.84,-0.84,0.16 \ -1.25,-0.25,0.75 \ -0.75,0.25,1.25 \ 0.16,0.84,1.84)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(4) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-1.84,-0.84,0.16 \ -1.25,-0.25,0.75 \ -0.75,0.25,1.25 \ 0.16,0.84,1.84)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(4) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_2A_50.csv", replace
* Scenario 2B : H_0 is FALSE / Effect size = 0.2
di "SCENARIO 2B - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-1.84,-0.84,0.16 \ -1.25,-0.25,0.75 \ -0.75,0.25,1.25 \ 0.16,0.84,1.84)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(4) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-1.84,-0.84,0.16 \ -1.25,-0.25,0.75 \ -0.75,0.25,1.25 \ 0.16,0.84,1.84)
qui simirt, nbobs(`Nn') mu(0.2) cov(1) dim(4) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_2B_50.csv", replace
* Scenario 2C : H_0 is FALSE / Effect size = 0.4
di "SCENARIO 2C - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-1.84,-0.84,0.16 \ -1.25,-0.25,0.75 \ -0.75,0.25,1.25 \ 0.16,0.84,1.84)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(4) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-1.84,-0.84,0.16 \ -1.25,-0.25,0.75 \ -0.75,0.25,1.25 \ 0.16,0.84,1.84)
qui simirt, nbobs(`Nn') mu(0.4) cov(1) dim(4) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_2C_50.csv", replace
* Scenario 2D : H_0 is FALSE / Effect size = -0.2
di "SCENARIO 2D - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-1.84,-0.84,0.16 \ -1.25,-0.25,0.75 \ -0.75,0.25,1.25 \ 0.16,0.84,1.84)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(4) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-1.84,-0.84,0.16 \ -1.25,-0.25,0.75 \ -0.75,0.25,1.25 \ 0.16,0.84,1.84)
qui simirt, nbobs(`Nn') mu(-0.2) cov(1) dim(4) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_2D_50.csv", replace
* Scenario 2E : H_0 is FALSE / Effect size = -0.4
di "SCENARIO 2E - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-1.84,-0.84,0.16 \ -1.25,-0.25,0.75 \ -0.75,0.25,1.25 \ 0.16,0.84,1.84)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(4) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-1.84,-0.84,0.16 \ -1.25,-0.25,0.75 \ -0.75,0.25,1.25 \ 0.16,0.84,1.84)
qui simirt, nbobs(`Nn') mu(-0.4) cov(1) dim(4) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_2E_50.csv", replace
**------------------------------------------------------------------------------------**
** Scenario 3: J = 7 items / M = 2 modalities / N=50 per group / TT=treatment variable
* Scenario 3A : H_0 is TRUE
di "SCENARIO 3A - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-1.15 \ -0.67 \ -0.32 \ 0 \ 0.32 \ 0.67 \ 1.15)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(7) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-1.15 \ -0.67 \ -0.32 \ 0 \ 0.32 \ 0.67 \ 1.15)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(7) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_3A_50.csv", replace
* Scenario 3B : H_0 is FALSE / Effect size = 0.2
di "SCENARIO 3B - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-1.15 \ -0.67 \ -0.32 \ 0 \ 0.32 \ 0.67 \ 1.15)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(7) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-1.15 \ -0.67 \ -0.32 \ 0 \ 0.32 \ 0.67 \ 1.15)
qui simirt, nbobs(`Nn') mu(0.2) cov(1) dim(7) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_3B_50.csv", replace
* Scenario 3C : H_0 is FALSE / Effect size = 0.4
di "SCENARIO 3C - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-1.15 \ -0.67 \ -0.32 \ 0 \ 0.32 \ 0.67 \ 1.15)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(7) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-1.15 \ -0.67 \ -0.32 \ 0 \ 0.32 \ 0.67 \ 1.15)
qui simirt, nbobs(`Nn') mu(0.4) cov(1) dim(7) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_3C_50.csv", replace
* Scenario 3D : H_0 is FALSE / Effect size = -0.2
di "SCENARIO 3D - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-1.15 \ -0.67 \ -0.32 \ 0 \ 0.32 \ 0.67 \ 1.15)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(7) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-1.15 \ -0.67 \ -0.32 \ 0 \ 0.32 \ 0.67 \ 1.15)
qui simirt, nbobs(`Nn') mu(-0.2) cov(1) dim(7) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_3D_50.csv", replace
* Scenario 3E : H_0 is FALSE / Effect size = -0.4
di "SCENARIO 3E - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-1.15 \ -0.67 \ -0.32 \ 0 \ 0.32 \ 0.67 \ 1.15)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(7) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-1.15 \ -0.67 \ -0.32 \ 0 \ 0.32 \ 0.67 \ 1.15)
qui simirt, nbobs(`Nn') mu(-0.4) cov(1) dim(7) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_3E_50.csv", replace
**------------------------------------------------------------------------------------**
** Scenario 4: J = 7 items / M = 4 modalities / N=50 per group / TT=treatment variable
* Scenario 4A : H_0 is TRUE
di "SCENARIO 4A - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-2.15,-1.15,-0.15 \ -1.67,-0.67,0.33 \ -1.32,-0.32,0.68 \ -1,0,1 \ -0.68,0.32,1.32 \ -0.33,0.67,1.67 \ 0.15,1.15,2.15)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(7) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-2.15,-1.15,-0.15 \ -1.67,-0.67,0.33 \ -1.32,-0.32,0.68 \ -1,0,1 \ -0.68,0.32,1.32 \ -0.33,0.67,1.67 \ 0.15,1.15,2.15)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(7) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_4A_50.csv", replace
* Scenario 4B : H_0 is FALSE / Effect size = 0.2
di "SCENARIO 4B - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-2.15,-1.15,-0.15 \ -1.67,-0.67,0.33 \ -1.32,-0.32,0.68 \ -1,0,1 \ -0.68,0.32,1.32 \ -0.33,0.67,1.67 \ 0.15,1.15,2.15)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(7) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-2.15,-1.15,-0.15 \ -1.67,-0.67,0.33 \ -1.32,-0.32,0.68 \ -1,0,1 \ -0.68,0.32,1.32 \ -0.33,0.67,1.67 \ 0.15,1.15,2.15)
qui simirt, nbobs(`Nn') mu(0.2) cov(1) dim(7) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_4B_50.csv", replace
* Scenario 4C : H_0 is FALSE / Effect size = 0.4
di "SCENARIO 4C - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-2.15,-1.15,-0.15 \ -1.67,-0.67,0.33 \ -1.32,-0.32,0.68 \ -1,0,1 \ -0.68,0.32,1.32 \ -0.33,0.67,1.67 \ 0.15,1.15,2.15)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(7) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-2.15,-1.15,-0.15 \ -1.67,-0.67,0.33 \ -1.32,-0.32,0.68 \ -1,0,1 \ -0.68,0.32,1.32 \ -0.33,0.67,1.67 \ 0.15,1.15,2.15)
qui simirt, nbobs(`Nn') mu(0.4) cov(1) dim(7) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_4C_50.csv", replace
* Scenario 4D : H_0 is FALSE / Effect size = -0.2
di "SCENARIO 4D - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-2.15,-1.15,-0.15 \ -1.67,-0.67,0.33 \ -1.32,-0.32,0.68 \ -1,0,1 \ -0.68,0.32,1.32 \ -0.33,0.67,1.67 \ 0.15,1.15,2.15)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(7) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-2.15,-1.15,-0.15 \ -1.67,-0.67,0.33 \ -1.32,-0.32,0.68 \ -1,0,1 \ -0.68,0.32,1.32 \ -0.33,0.67,1.67 \ 0.15,1.15,2.15)
qui simirt, nbobs(`Nn') mu(-0.2) cov(1) dim(7) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_4D_50.csv", replace
* Scenario 4E : H_0 is FALSE / Effect size = -0.4
di "SCENARIO 4E - N=50"
forvalues replication = 1/1000 {
if mod(`replication',100)==0 {
di "replication = `replication'"
}
mat D= (-2.15,-1.15,-0.15 \ -1.67,-0.67,0.33 \ -1.32,-0.32,0.68 \ -1,0,1 \ -0.68,0.32,1.32 \ -0.33,0.67,1.67 \ 0.15,1.15,2.15)
qui simirt, nbobs(`Nn') mu(0) cov(1) dim(7) pcm(D) clear
qui gen TT = 0
tempfile grp0
qui save `grp0',replace
mat D= (-2.15,-1.15,-0.15 \ -1.67,-0.67,0.33 \ -1.32,-0.32,0.68 \ -1,0,1 \ -0.68,0.32,1.32 \ -0.33,0.67,1.67 \ 0.15,1.15,2.15)
qui simirt, nbobs(`Nn') mu(-0.4) cov(1) dim(7) pcm(D) clear
qui gen TT = 1
tempfile grp1
qui save `grp1',replace
clear
use `grp0'
qui append using `grp1'
drop id
qui gen id = _n
order(id)
qui gen replication = `replication'
if `replication'==1{
tempfile data
qui save `data'
}
else{
qui append using `data'
qui save `data',replace
}
}
export delimited using "`path'/scenario_4E_50.csv", replace