File restructure + reproduction guide
This commit is contained in:
438
RProject/Scripts/functions/desc_analysis_dif.R
Normal file
438
RProject/Scripts/functions/desc_analysis_dif.R
Normal file
@ -0,0 +1,438 @@
|
||||
##############################################################################
|
||||
#----------------------------------------------------------------------------#
|
||||
############################ BOXPLOTS H0 SCENARIOS ###########################
|
||||
#----------------------------------------------------------------------------#
|
||||
##############################################################################
|
||||
|
||||
## Proportion of rejected h0 per dif value in h0 scenarios (A) by DIF size
|
||||
|
||||
res.null <- res.dat[res.dat$eff.size==0,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,1),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
|
||||
res.null0 <- res.dat[res.dat$eff.size==0 & res.dat$dif.size==0,]
|
||||
points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='gray',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size==0 & res.dat$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size==0 & res.dat$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
## Proportion of rejected h0 per dif value in h0 scenarios (A) by DIF size
|
||||
|
||||
par(mfrow=c(2,2))
|
||||
|
||||
# 1 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 2 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 3 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
par(mfrow=c(1,1))
|
||||
|
||||
## Proportion of rejected h0 per dif value in h0 scenarios (A) by DIF size (1 item) J=4
|
||||
|
||||
par(mfrow=c(2,2))
|
||||
|
||||
# 1 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1 & res.dat$J==4,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 2 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$J==4,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
par(mfrow=c(1,1))
|
||||
|
||||
|
||||
## Proportion of rejected h0 per dif value in h0 scenarios (A) by DIF size (1 item) J=7
|
||||
|
||||
par(mfrow=c(2,2))
|
||||
|
||||
# 2 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$J==7,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 3 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3 & res.dat$J==7,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
par(mfrow=c(1,1))
|
||||
|
||||
############# By N
|
||||
|
||||
####### N=100
|
||||
|
||||
|
||||
par(mfrow=c(2,2))
|
||||
|
||||
# 1 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1 & res.dat$N==100,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 2 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$N==100,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 3 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3 & res.dat$N==100,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
par(mfrow=c(1,1))
|
||||
|
||||
|
||||
####### N=200
|
||||
|
||||
|
||||
par(mfrow=c(2,2))
|
||||
|
||||
# 1 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1 & res.dat$N==200,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 2 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$N==200,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 3 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3 & res.dat$N==200,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
par(mfrow=c(1,1))
|
||||
|
||||
|
||||
####### N=300
|
||||
|
||||
|
||||
par(mfrow=c(2,2))
|
||||
|
||||
# 1 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1 & res.dat$N==300,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 2 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$N==300,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 3 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3 & res.dat$N==300,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
par(mfrow=c(1,1))
|
||||
|
||||
##############################################################################
|
||||
#----------------------------------------------------------------------------#
|
||||
############################ BOXPLOTS H1 SCENARIOS ###########################
|
||||
#----------------------------------------------------------------------------#
|
||||
##############################################################################
|
||||
|
||||
#### CALCULER LA PUISSANCE THEORIQUE AVEC RASCHPOWER
|
||||
|
||||
## Proportion of rejected h0 per dif value in h1 scenarios by DIF size // eff.size positive
|
||||
|
||||
res.null <- res.dat[res.dat$eff.size>0,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
|
||||
res.null0 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0,]
|
||||
points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(5,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
|
||||
############# By N // EFF SIZE POSITIVE
|
||||
|
||||
####### N=100
|
||||
|
||||
res.null <- res.dat[res.dat$eff.size>0 & res.dat$N==100,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
|
||||
res.null0 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0 & res.dat$N==100,]
|
||||
points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.3 & res.dat$N==100,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.5 & res.dat$N==100,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.3 & res.dat$N==100,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.5 & res.dat$N==100,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(5,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
####### N=300 // DIF à 0.5 - QUELS SONT LES SCENARIOS EN HAUT
|
||||
|
||||
res.null <- res.dat[res.dat$eff.size>0 & res.dat$N==300,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
|
||||
res.null0 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0 & res.dat$N==300,]
|
||||
points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.3 & res.dat$N==300,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.5 & res.dat$N==300,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.3 & res.dat$N==300,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.5 & res.dat$N==300,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(5,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
|
||||
############# By N // EFF SIZE NEGATIVE
|
||||
|
||||
####### N=100
|
||||
|
||||
res.null <- res.dat[res.dat$eff.size<0 & res.dat$N==100,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
|
||||
res.null0 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==0 & res.dat$N==100,]
|
||||
points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==-0.3 & res.dat$N==100,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==-0.5 & res.dat$N==100,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
|
||||
####### N=300 // DIF à 0.5 - QUELS SONT LES SCENARIOS EN HAUT
|
||||
|
||||
res.null <- res.dat[res.dat$eff.size<0 & res.dat$N==300,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
|
||||
res.null0 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==0 & res.dat$N==300,]
|
||||
points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==-0.3 & res.dat$N==300,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==-0.5 & res.dat$N==300,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
|
||||
##############################################################################
|
||||
#----------------------------------------------------------------------------#
|
||||
########################## SYSTEMATIC ERROR BOXPLOTS #########################
|
||||
#----------------------------------------------------------------------------#
|
||||
##############################################################################
|
||||
|
||||
# Overall
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat,col=c(2,3),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat[res.dat$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
|
||||
|
||||
# J=4
|
||||
res.dat.temp <- res.dat[res.dat$J==4,]
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
|
||||
|
||||
|
||||
# J=7
|
||||
res.dat.temp <- res.dat[res.dat$J==7,]
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
|
||||
|
||||
# J=4 / 1 DIF
|
||||
res.dat.temp <- res.dat[res.dat$J==4 & res.dat$nb.dif==1,]
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
|
||||
|
||||
# J=4 / 2 DIF
|
||||
res.dat.temp <- res.dat[res.dat$J==4 & res.dat$nb.dif==2,]
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
|
||||
|
||||
|
||||
# J=7 / 2 DIF
|
||||
|
||||
res.dat.temp <- res.dat[res.dat$J==7 & res.dat$nb.dif==2,]
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
|
||||
|
||||
|
||||
# J=7 / 3 DIF
|
||||
|
||||
res.dat.temp <- res.dat[res.dat$J==7 & res.dat$nb.dif==3,]
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
|
||||
|
||||
|
||||
##############################################################################
|
||||
#----------------------------------------------------------------------------#
|
||||
########################## BETA SIGN CHANGE BOXPLOTS #########################
|
||||
#----------------------------------------------------------------------------#
|
||||
##############################################################################
|
||||
|
||||
# Overall
|
||||
|
||||
boxplot(beta.same.sign.truebeta.p~dif.size,data=res.dat,col=c(2,3),xlab='DIF size',
|
||||
ylab='Proportion of estimates with same sign as true value in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat[res.dat$dif.size==-0.5,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0.5,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==-0.3,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0.3,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
|
||||
|
||||
# Overall // H0 rejected
|
||||
|
||||
boxplot(beta.same.sign.truebeta.signif.p~dif.size,data=res.dat,col=c(2,3),xlab='DIF size',
|
||||
ylab='Proportion of estimates with same sign as true value in target scenario',main='When H0 rejected',ylim=c(0,1))
|
||||
res.null3 <- res.dat[res.dat$dif.size==-0.5,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0.5,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==-0.3,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0.3,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
|
||||
475
RProject/Scripts/functions/desc_analysis_nodif.R
Normal file
475
RProject/Scripts/functions/desc_analysis_nodif.R
Normal file
@ -0,0 +1,475 @@
|
||||
##############################################################################
|
||||
#----------------------------------------------------------------------------#
|
||||
############################ BOXPLOTS H0 SCENARIOS ###########################
|
||||
#----------------------------------------------------------------------------#
|
||||
##############################################################################
|
||||
|
||||
## Proportion of rejected h0 per dif value in h0 scenarios (A) by DIF size
|
||||
|
||||
res.null <- res.dat[res.dat$eff.size==0,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,1),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
|
||||
res.null0 <- res.dat[res.dat$eff.size==0 & res.dat$dif.size==0,]
|
||||
points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='gray',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size==0 & res.dat$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size==0 & res.dat$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
## Proportion of rejected h0 per dif value in h0 scenarios (A) by DIF size
|
||||
|
||||
par(mfrow=c(2,2))
|
||||
|
||||
# 0 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==0,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='No DIF',ylim=c(0,1))
|
||||
points(y=res.null$h0.rejected.p,x=rep(1,nrow(res.null)),col='#590b0c',pch=3)
|
||||
|
||||
# 1 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 2 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 3 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
par(mfrow=c(1,1))
|
||||
|
||||
## Proportion of rejected h0 per dif value in h0 scenarios (A) by DIF size (1 item) J=4
|
||||
|
||||
par(mfrow=c(2,2))
|
||||
|
||||
# 0 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==0 & res.dat$J==4,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='No DIF',ylim=c(0,1))
|
||||
points(y=res.null$h0.rejected.p,x=rep(1,nrow(res.null)),col='#590b0c',pch=3)
|
||||
|
||||
# 1 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1 & res.dat$J==4,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 2 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$J==4,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
par(mfrow=c(1,1))
|
||||
|
||||
|
||||
## Proportion of rejected h0 per dif value in h0 scenarios (A) by DIF size (1 item) J=7
|
||||
|
||||
par(mfrow=c(2,2))
|
||||
|
||||
# 0 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==0 & res.dat$J==7,]
|
||||
nrow(res.null)
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='No DIF',ylim=c(0,1))
|
||||
points(y=res.null$h0.rejected.p,x=rep(1,nrow(res.null)),col='#590b0c',pch=3)
|
||||
|
||||
# 2 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$J==7,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 3 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3 & res.dat$J==7,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
par(mfrow=c(1,1))
|
||||
|
||||
############# By N
|
||||
|
||||
####### N=100
|
||||
|
||||
|
||||
par(mfrow=c(2,2))
|
||||
|
||||
# 0 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==0 & res.dat$N==100,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='No DIF',ylim=c(0,1))
|
||||
points(y=res.null$h0.rejected.p,x=rep(1,nrow(res.null)),col='#590b0c',pch=3)
|
||||
|
||||
# 1 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1 & res.dat$N==100,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 2 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$N==100,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 3 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3 & res.dat$N==100,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
par(mfrow=c(1,1))
|
||||
|
||||
|
||||
####### N=200
|
||||
|
||||
|
||||
par(mfrow=c(2,2))
|
||||
|
||||
# 0 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==0 & res.dat$N==200,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='No DIF',ylim=c(0,1))
|
||||
points(y=res.null$h0.rejected.p,x=rep(1,nrow(res.null)),col='#590b0c',pch=3)
|
||||
|
||||
# 1 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1 & res.dat$N==200,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 2 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$N==200,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 3 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3 & res.dat$N==200,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
par(mfrow=c(1,1))
|
||||
|
||||
|
||||
####### N=300
|
||||
|
||||
|
||||
par(mfrow=c(2,2))
|
||||
|
||||
# 0 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==0 & res.dat$N==300,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='No DIF',ylim=c(0,1))
|
||||
points(y=res.null$h0.rejected.p,x=rep(1,nrow(res.null)),col='#590b0c',pch=3)
|
||||
|
||||
# 1 item
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==1 & res.dat$N==300,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 1 item',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 2 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==2 & res.dat$N==300,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 2 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
# 3 items
|
||||
res.null <- res.dat[res.dat$eff.size==0 & res.dat$nb.dif==3 & res.dat$N==300,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2),xlab='DIF size',
|
||||
ylab='H0 rejection proportion in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.null[res.null$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.null[res.null$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
par(mfrow=c(1,1))
|
||||
|
||||
##############################################################################
|
||||
#----------------------------------------------------------------------------#
|
||||
############################ BOXPLOTS H1 SCENARIOS ###########################
|
||||
#----------------------------------------------------------------------------#
|
||||
##############################################################################
|
||||
|
||||
#### CALCULER LA PUISSANCE THEORIQUE AVEC RASCHPOWER
|
||||
|
||||
## Proportion of rejected h0 per dif value in h1 scenarios by DIF size // eff.size positive
|
||||
|
||||
res.null <- res.dat[res.dat$eff.size>0,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
|
||||
res.null0 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0,]
|
||||
points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.3,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.5,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(5,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
|
||||
############# By N // EFF SIZE POSITIVE
|
||||
|
||||
####### N=100
|
||||
|
||||
res.null <- res.dat[res.dat$eff.size>0 & res.dat$N==100,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
|
||||
res.null0 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0 & res.dat$N==100,]
|
||||
points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.3 & res.dat$N==100,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.5 & res.dat$N==100,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.3 & res.dat$N==100,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.5 & res.dat$N==100,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(5,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
####### N=300 // DIF à 0.5 - QUELS SONT LES SCENARIOS EN HAUT
|
||||
|
||||
res.null <- res.dat[res.dat$eff.size>0 & res.dat$N==300,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
|
||||
res.null0 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0 & res.dat$N==300,]
|
||||
points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.3 & res.dat$N==300,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==-0.5 & res.dat$N==300,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.3 & res.dat$N==300,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size>0 & res.dat$dif.size==0.5 & res.dat$N==300,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(5,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
|
||||
############# By N // EFF SIZE NEGATIVE
|
||||
|
||||
####### N=100
|
||||
|
||||
res.null <- res.dat[res.dat$eff.size<0 & res.dat$N==100,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
|
||||
res.null0 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==0 & res.dat$N==100,]
|
||||
points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==-0.3 & res.dat$N==100,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==-0.5 & res.dat$N==100,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
|
||||
####### N=300 // DIF à 0.5 - QUELS SONT LES SCENARIOS EN HAUT
|
||||
|
||||
res.null <- res.dat[res.dat$eff.size<0 & res.dat$N==300,]
|
||||
boxplot(h0.rejected.p~dif.size,data=res.null,col=c(3,2,4,2,3),xlab='DIF size',ylab='H0 rejection proportion in target scenario')
|
||||
res.null0 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==0 & res.dat$N==300,]
|
||||
points(y=res.null0$h0.rejected.p,x=rep(3,nrow(res.null0)),col='darkblue',pch=3)
|
||||
res.null3 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==-0.3 & res.dat$N==300,]
|
||||
points(y=res.null3$h0.rejected.p,x=rep(2,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null5 <- res.dat[res.dat$eff.size<0 & res.dat$dif.size==-0.5 & res.dat$N==300,]
|
||||
points(y=res.null5$h0.rejected.p,x=rep(1,nrow(res.null5)),col='#053305',pch=3)
|
||||
|
||||
|
||||
##############################################################################
|
||||
#----------------------------------------------------------------------------#
|
||||
########################## SYSTEMATIC ERROR BOXPLOTS #########################
|
||||
#----------------------------------------------------------------------------#
|
||||
##############################################################################
|
||||
|
||||
# Overall
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat,col=c(2,3),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat[res.dat$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
|
||||
|
||||
# J=4
|
||||
res.dat.temp <- res.dat[res.dat$J==4,]
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
|
||||
|
||||
|
||||
# J=7
|
||||
res.dat.temp <- res.dat[res.dat$J==7,]
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
|
||||
|
||||
# J=4 / 1 DIF
|
||||
res.dat.temp <- res.dat[res.dat$J==4 & res.dat$nb.dif==1,]
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
|
||||
|
||||
# J=4 / 2 DIF
|
||||
res.dat.temp <- res.dat[res.dat$J==4 & res.dat$nb.dif==2,]
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
|
||||
|
||||
|
||||
# J=7 / 2 DIF
|
||||
|
||||
res.dat.temp <- res.dat[res.dat$J==7 & res.dat$nb.dif==2,]
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
|
||||
|
||||
|
||||
# J=7 / 3 DIF
|
||||
|
||||
res.dat.temp <- res.dat[res.dat$J==7 & res.dat$nb.dif==3,]
|
||||
|
||||
boxplot(true.value.in.ci.p~dif.size,data=res.dat.temp,col=c(2,3,3,2),xlab='DIF size',
|
||||
ylab='Proportion of true beta value in CI in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.5,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(4,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==-0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat.temp[res.dat.temp$dif.size==0.3,]
|
||||
points(y=res.null3$true.value.in.ci.p,x=rep(3,nrow(res.null3)),col='#053305',pch=3)
|
||||
|
||||
|
||||
##############################################################################
|
||||
#----------------------------------------------------------------------------#
|
||||
########################## BETA SIGN CHANGE BOXPLOTS #########################
|
||||
#----------------------------------------------------------------------------#
|
||||
##############################################################################
|
||||
|
||||
# Overall
|
||||
|
||||
boxplot(beta.same.sign.truebeta.p~dif.size,data=res.dat,col=c(2,3),xlab='DIF size',
|
||||
ylab='Proportion of estimates with same sign as true value in target scenario',main='DIF on 3 items',ylim=c(0,1))
|
||||
res.null3 <- res.dat[res.dat$dif.size==-0.5,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0.5,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==-0.3,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0.3,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
|
||||
|
||||
# Overall // H0 rejected
|
||||
|
||||
boxplot(beta.same.sign.truebeta.signif.p~dif.size,data=res.dat,col=c(2,3),xlab='DIF size',
|
||||
ylab='Proportion of estimates with same sign as true value in target scenario',main='When H0 rejected',ylim=c(0,1))
|
||||
res.null3 <- res.dat[res.dat$dif.size==-0.5,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(1,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0.5,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(5,nrow(res.null3)),col='#590b0c',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==-0.3,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(2,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0.3,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(4,nrow(res.null3)),col='#053305',pch=3)
|
||||
res.null3 <- res.dat[res.dat$dif.size==0,]
|
||||
points(y=res.null3$beta.same.sign.truebeta.signif.p,x=rep(3,nrow(res.null3)),col='gray',pch=3)
|
||||
293
RProject/Scripts/functions/power_analysis.R
Normal file
293
RProject/Scripts/functions/power_analysis.R
Normal 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')
|
||||
270
RProject/Scripts/functions/resali.R
Normal file
270
RProject/Scripts/functions/resali.R
Normal file
@ -0,0 +1,270 @@
|
||||
library(TAM)
|
||||
|
||||
resali <- function(df=NULL,items=NULL,group=NULL,verbose=T) {
|
||||
if (verbose) {
|
||||
cat('-----------------------------------------------------------\n')
|
||||
cat('COMPUTING INITIAL PCM\n')
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
nbitems <- length(items)
|
||||
nbitems_o <- nbitems
|
||||
items_n <- paste0('item',items)
|
||||
resp <- df[,items_n]
|
||||
grp <- df[,group]
|
||||
pcm_initial <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
|
||||
dat <- df
|
||||
dat$score <- rowSums(dat[,items_n])
|
||||
nqt <- ifelse(length(unique(quantile(dat$score,seq(0,1,0.2))))==6,5,length(unique(quantile(dat$score,seq(0,1,0.2))))-1)
|
||||
while (length(unique(quantile(dat$score,seq(0,1,1/nqt))))!=nqt+1) {
|
||||
nqt <- nqt-1
|
||||
}
|
||||
dat$score_q5 <- cut(dat$score,unique(quantile(dat$score,seq(0,1,1/nqt))),labels=1:nqt,include.lowest=T)
|
||||
res.anova <- rep(NA,nbitems)
|
||||
pval <- rep(NA,nbitems)
|
||||
fval <- rep(NA,nbitems)
|
||||
for (i in items) {
|
||||
dat[,paste0('res_',i)] <- IRT.residuals(pcm_initial)$stand_residuals[,i]
|
||||
res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~TT*score_q5,data=dat))
|
||||
pval[i] <- res.anova[[i]][1,"Pr(>F)"]
|
||||
fval[i] <- res.anova[[i]][1,'F value']
|
||||
}
|
||||
if (verbose) {
|
||||
cat('DONE\n')
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
res.items <- c()
|
||||
res.uniform <- c()
|
||||
k <- 1
|
||||
while(any(pval<0.05/(nbitems_o*3))) {
|
||||
k <- k+1
|
||||
if (verbose) {
|
||||
cat(paste('COMPUTING STEP',k,'\n'))
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
res.item <- gsub("[a-z]", "",colnames(resp)[which.max(fval)])
|
||||
res.items <- c(res.items,res.item)
|
||||
res.uni <- res.anova[[which.max(fval)]][3,"Pr(>F)"]>0.05
|
||||
res.uniform <- c(res.uniform,res.uni)
|
||||
items_n <- c(items_n[items_n!=paste0('item',res.item)],paste0("item",res.item,c("noTT","TT")))
|
||||
dat[,paste0("item",res.item,'TT')] <- dat[dat$TT==1,paste0('item',res.item)]
|
||||
dat[,paste0("item",res.item,'noTT')] <- dat[dat$TT==0,paste0('item',res.item)]
|
||||
resp <- dat[,items_n]
|
||||
grp <- dat[,group]
|
||||
pcm_while <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
|
||||
nbitems <- length(items_n)
|
||||
res.anova <- rep(NA,nbitems)
|
||||
pval <- rep(NA,nbitems)
|
||||
fval <- rep(NA,nbitems)
|
||||
for (i in 1:nbitems) {
|
||||
dat[,paste0('res_',i)] <- IRT.residuals(pcm_while)$stand_residuals[,i]
|
||||
res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~TT*score_q5,data=dat))
|
||||
pval[i] <- res.anova[[i]][1,"Pr(>F)"]
|
||||
fval[i] <- res.anova[[i]][1,'F value']
|
||||
}
|
||||
if (verbose) {
|
||||
cat('DONE\n')
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
}
|
||||
if (verbose) {
|
||||
cat("DETECTED DIF ITEMS\n")
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
if (length(res.items>0)) {
|
||||
results <- data.frame(dif.items=res.items,
|
||||
uniform=1*res.uniform)
|
||||
return(results)
|
||||
}
|
||||
else {
|
||||
if (verbose) {
|
||||
cat("No DIF was detected\n")
|
||||
}
|
||||
return(NULL)
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
resali_BF2 <- function(df=NULL,items=NULL,group=NULL,verbose=T) {
|
||||
if (verbose) {
|
||||
cat('-----------------------------------------------------------\n')
|
||||
cat('COMPUTING INITIAL PCM\n')
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
nbitems <- length(items)
|
||||
nbitems_o <- nbitems
|
||||
items_n <- paste0('item',items)
|
||||
resp <- df[,items_n]
|
||||
grp <- df[,group]
|
||||
pcm_initial <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
|
||||
dat <- df
|
||||
dat$score <- rowSums(dat[,items_n])
|
||||
nqt <- ifelse(length(unique(quantile(dat$score,seq(0,1,0.2))))==6,5,length(unique(quantile(dat$score,seq(0,1,0.2))))-1)
|
||||
dat$score_q5 <- cut(dat$score,unique(quantile(dat$score,seq(0,1,1/nqt))),labels=1:nqt,include.lowest=T)
|
||||
res.anova <- rep(NA,nbitems)
|
||||
pval <- rep(NA,nbitems)
|
||||
fval <- rep(NA,nbitems)
|
||||
for (i in items) {
|
||||
dat[,paste0('res_',i)] <- rowMeans(predict(pcm_initial)$stand.resid[,,i])
|
||||
res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~TT*score_q5,data=dat))
|
||||
pval[i] <- res.anova[[i]][1,"Pr(>F)"]
|
||||
fval[i] <- res.anova[[i]][1,'F value']
|
||||
}
|
||||
if (verbose) {
|
||||
cat('DONE\n')
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
res.items <- c()
|
||||
res.uniform <- c()
|
||||
k <- 1
|
||||
while(any(pval<0.05/(nbitems_o-k+1))) {
|
||||
k <- k+1
|
||||
if (verbose) {
|
||||
cat(paste('COMPUTING STEP',k,'\n'))
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
res.item <- gsub("[a-z]", "",colnames(resp)[which.max(fval)])
|
||||
res.items <- c(res.items,res.item)
|
||||
res.uni <- res.anova[[which.max(fval)]][3,"Pr(>F)"]>0.05
|
||||
res.uniform <- c(res.uniform,res.uni)
|
||||
items_n <- c(items_n[items_n!=paste0('item',res.item)],paste0("item",res.item,c("noTT","TT")))
|
||||
dat[,paste0("item",res.item,'TT')] <- dat[dat$TT==1,paste0('item',res.item)]
|
||||
dat[,paste0("item",res.item,'noTT')] <- dat[dat$TT==0,paste0('item',res.item)]
|
||||
resp <- dat[,items_n]
|
||||
grp <- dat[,group]
|
||||
pcm_while <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
|
||||
nbitems <- length(items_n)
|
||||
res.anova <- rep(NA,nbitems)
|
||||
pval <- rep(NA,nbitems)
|
||||
fval <- rep(NA,nbitems)
|
||||
for (i in 1:nbitems) {
|
||||
dat[,paste0('res_',i)] <- rowMeans(predict(pcm_while)$stand.resid[,,i])
|
||||
res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~TT*score_q5,data=dat))
|
||||
pval[i] <- res.anova[[i]][1,"Pr(>F)"]
|
||||
fval[i] <- res.anova[[i]][1,'F value']
|
||||
}
|
||||
if (verbose) {
|
||||
cat('DONE\n')
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
}
|
||||
if (verbose) {
|
||||
cat("DETECTED DIF ITEMS\n")
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
if (length(res.items>0)) {
|
||||
results <- data.frame(dif.items=res.items,
|
||||
uniform=1*res.uniform)
|
||||
return(results)
|
||||
}
|
||||
else {
|
||||
if (verbose) {
|
||||
cat("No DIF was detected\n")
|
||||
}
|
||||
return(NULL)
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
resali_LRT <- function(df=NULL,items=NULL,group=NULL,verbose=T) {
|
||||
if (verbose) {
|
||||
cat('-----------------------------------------------------------\n')
|
||||
cat('COMPUTING INITIAL PCM\n')
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
nbitems <- length(items)
|
||||
nbitems_o <- nbitems
|
||||
items_n <- paste0('item',items)
|
||||
resp <- df[,items_n]
|
||||
grp <- df[,group]
|
||||
items_n_alt <- paste0(items_n,c("_TT","_noTT"))
|
||||
for (i in items_n) {
|
||||
df[df$TT==0,paste0(i,"_noTT")] <- df[df$TT==0,i]
|
||||
df[df$TT==1,paste0(i,"_TT")] <- df[df$TT==1,i]
|
||||
}
|
||||
resp_alt <- df[,items_n_alt]
|
||||
pcm_initial <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
|
||||
pcm_alt <- TAM::tam.mml(resp=resp_alt,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
|
||||
dat <- df
|
||||
dat$score <- rowSums(dat[,items_n])
|
||||
nqt <- ifelse(length(unique(quantile(dat$score,seq(0,1,0.2))))==6,5,length(unique(quantile(dat$score,seq(0,1,0.2))))-1)
|
||||
dat$score_q5 <- cut(dat$score,unique(quantile(dat$score,seq(0,1,1/nqt))),labels=1:nqt,include.lowest=T)
|
||||
res.anova <- rep(NA,nbitems)
|
||||
pval <- rep(NA,nbitems)
|
||||
fval <- rep(NA,nbitems)
|
||||
for (i in items) {
|
||||
dat[,paste0('res_',i)] <- rowMeans(predict(pcm_initial)$stand.resid[,,i])
|
||||
res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~TT*score_q5,data=dat))
|
||||
pval[i] <- res.anova[[i]][1,"Pr(>F)"]
|
||||
fval[i] <- res.anova[[i]][1,'F value']
|
||||
}
|
||||
if (verbose) {
|
||||
cat('DONE\n')
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
res.items <- c()
|
||||
res.uniform <- c()
|
||||
k <- 1
|
||||
if (anova(pcm_initial,pcm_alt)$p[1]<0.05) {
|
||||
while(any(pval<0.05/nbitems_o)) {
|
||||
k <- k+1
|
||||
if (verbose) {
|
||||
cat(paste('COMPUTING STEP',k,'\n'))
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
res.item <- gsub("[a-z]", "",colnames(resp)[which.max(fval)])
|
||||
res.items <- c(res.items,res.item)
|
||||
res.uni <- res.anova[[which.max(fval)]][3,"Pr(>F)"]>0.05
|
||||
res.uniform <- c(res.uniform,res.uni)
|
||||
items_n <- c(items_n[items_n!=paste0('item',res.item)],paste0("item",res.item,c("noTT","TT")))
|
||||
dat[,paste0("item",res.item,'TT')] <- dat[dat$TT==1,paste0('item',res.item)]
|
||||
dat[,paste0("item",res.item,'noTT')] <- dat[dat$TT==0,paste0('item',res.item)]
|
||||
resp <- dat[,items_n]
|
||||
grp <- dat[,group]
|
||||
pcm_while <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = "PCM",est.variance = T,verbose=F)
|
||||
nbitems <- length(items_n)
|
||||
res.anova <- rep(NA,nbitems)
|
||||
pval <- rep(NA,nbitems)
|
||||
fval <- rep(NA,nbitems)
|
||||
for (i in 1:nbitems) {
|
||||
dat[,paste0('res_',i)] <- rowMeans(predict(pcm_while)$stand.resid[,,i])
|
||||
res.anova[i] <- summary(aov(dat[,paste0('res_',i)]~TT*score_q5,data=dat))
|
||||
pval[i] <- res.anova[[i]][1,"Pr(>F)"]
|
||||
fval[i] <- res.anova[[i]][1,'F value']
|
||||
}
|
||||
if (verbose) {
|
||||
cat('DONE\n')
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
}
|
||||
if (verbose) {
|
||||
cat("DETECTED DIF ITEMS\n")
|
||||
cat('-----------------------------------------------------------\n')
|
||||
}
|
||||
if (length(res.items>0)) {
|
||||
results <- data.frame(dif.items=res.items,
|
||||
uniform=1*res.uniform)
|
||||
return(results)
|
||||
}
|
||||
else {
|
||||
if (verbose) {
|
||||
cat("No DIF was detected\n")
|
||||
}
|
||||
return(NULL)
|
||||
}
|
||||
}
|
||||
else {
|
||||
if (verbose) {
|
||||
cat("No DIF was detected\n")
|
||||
}
|
||||
return(NULL)
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user