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

Corrected bug in residuals analysis

Browse Source
main
Corentin Choisy 8 months ago
parent d65132a190
commit 7d9f634113
2 changed files with 128 additions and 14 deletions
Show Stats Download Patch File Download Diff File

134
RProject/Scripts/aggregation.R
Unescape Escape View File

@ -629,24 +629,39 @@ results <- c(results,results2)
compile_simulation2_resali <- function(scenario) {
name <- as.numeric(gsub("[^0-9.-]", "", substr(scenario,start=0,stop=2)))
if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50" & name>0) {
s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/resali-DIF/N50/',scenario,'_original.xls'))
s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Results/N50/',scenario,'_original.xls'))
}
if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name>0) {
s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/resali-DIF/N100/',scenario,'_original.xls'))
s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Results/N100/',scenario,'_original.xls'))
}
if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name>0) {
s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/resali-DIF/N200/',scenario,'_original.xls'))
s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Results/N200/',scenario,'_original.xls'))
}
if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name>0) {
s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/resali-DIF/N300/',scenario,'_original.xls'))
s <- read_excel(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/RESALI/Results/N300/',scenario,'_original.xls'))
}
J <- max(which(sapply(1:7,function(x) paste0('item',x) %in% colnames(s) | paste0('item',x,'_1') %in% colnames(s))))
M <- 1+sum(sapply(1:3,function(x) paste0('item1_',x) %in% colnames(s) ))
if (M==1) {M <- 2}
nb.dif <- max(which(sapply(1:3,function(x) paste0('dif',x) %in% colnames(s) | paste0('dif',x,'_1') %in% colnames(s))))
nb.dif.true <- ifelse(name<=4,0,ifelse(name<=8,1,ifelse(name<=16,2,3)))
if (name %in% c(3,4,13:20)) {
m.nb.dif.detect <- mean(ifelse(is.na(s$dif_1_1),0,
ifelse(is.na(s$dif_2_1),1,
ifelse(is.na(s$dif_3_1),2,
ifelse(is.na(s$dif_4_1),3,
ifelse(is.na(s$dif_5_1),4,
ifelse(is.na(s$dif_6_1),5,
ifelse(is.na(s$dif_7_1),6,7))))))))
}
if (!(name %in% c(3,4,13:20))) {
m.nb.dif.detect <- mean(ifelse(is.na(s$dif_1_1),0,
ifelse(is.na(s$dif_2_1),1,
ifelse(is.na(s$dif_3_1),2,
ifelse(is.na(s$dif_4_1),3,4)))))
}
if (J==4) {
if (M==2) {
a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4))
a <- data.frame(m.item1=mean(s$item1_1),m.item2=mean(s$item2_1),m.item3=mean(s$item3_1),m.item4=mean(s$item4_1))
} else {
a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
@ -656,8 +671,8 @@ compile_simulation2_resali <- function(scenario) {
}
} else {
if (M==2) {
a <- data.frame(m.item1=mean(s$item1),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4),
m.item5=mean(s$item5),m.item6=mean(s$item6),m.item7=mean(s$item7))
a <- data.frame(m.item1=mean(s$item1_1),m.item2=mean(s$item2_1),m.item3=mean(s$item3_1),m.item4=mean(s$item4_1),
m.item5=mean(s$item5_1),m.item6=mean(s$item6_1),m.item7=mean(s$item7_1))
} else {
a <- data.frame(m.item1_1=mean(s$item1_1),m.item1_2=mean(s$item1_2),m.item1_3=mean(s$item1_3),
m.item2_1=mean(s$item2_1),m.item2_2=mean(s$item2_2),m.item2_3=mean(s$item2_3),
@ -679,12 +694,106 @@ compile_simulation2_resali <- function(scenario) {
J=J,
M=M,
eff.size=eff.size,
nb.dif=nb.dif,
nb.dif=nb.dif.true,
m.nb.dif.detect=m.nb.dif.detect,
dif.size=dif.size
)
true.value.in.ci <- eff.size <= s$beta+1.96*s$se_beta & eff.size >= s$beta-1.96*s$se_beta
beta.same.sign.truebeta.p <- ifelse(rep(eff.size,nrow(s))==0,NA,(rep(eff.size,nrow(s))/s$beta)>0)
num.reject <- which((s$beta-1.96*s$se_beta)>0 | (s$beta+1.96*s$se_beta)<0)
dif.d <- mean(sapply(1:1000,function(x) any(!is.na(s[x,paste0("dif_",1:unique(b$J),"_1")]))))
if (nb.dif.true==0) {
prop.perfect <- NA
flexible.detect <- NA
moreflexible.detect <- NA
}
if (nb.dif.true==1 & unique(b$J)==4 & unique(b$M)==4) {
perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))==1,s[x,"dif_detect_unif_1"]==1 & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1")])
,0) )
prop.perfect <- mean(perfect.detection)
flexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))==1,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1")])
,0) )
flexible.detect <- mean(flexible.detect)
moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:4)]),0) )
moreflexible.detect <- mean(moreflexible.detect)
}
if (nb.dif.true==2 & unique(b$J)==4 & unique(b$M)==4) {
perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))==2,s[x,"dif_detect_unif_1"]==1 & s[x,"dif_detect_unif_2"]==1 & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2")])
,0) )
prop.perfect <- mean(perfect.detection)
flexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))==2,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2")])
,0) )
flexible.detect <- mean(flexible.detect)
moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:4)]) &
s[x,"real_dif_2"]%in%c(s[x,paste0("dif_detect_",1:4)]),0) )
moreflexible.detect <- mean(moreflexible.detect)
}
if (nb.dif.true==2 & unique(b$J)==7 & unique(b$M)==4) {
perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))==2,s[x,"dif_detect_unif_1"]==1 & s[x,"dif_detect_unif_2"]==1 & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2")])
,0) )
prop.perfect <- mean(perfect.detection)
flexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))==2,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2")])
,0) )
flexible.detect <- mean(flexible.detect)
moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:7)]) &
s[x,"real_dif_2"]%in%c(s[x,paste0("dif_detect_",1:7)]),0) )
moreflexible.detect <- mean(moreflexible.detect)
}
if (nb.dif.true==3 & unique(b$J)==7 & unique(b$M)==4) {
perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))==3,s[x,"dif_detect_unif_1"]==1 & s[x,"dif_detect_unif_2"]==1 & s[x,"dif_detect_unif_3"]==1 & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[3])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")])
,0) )
prop.perfect <- mean(perfect.detection)
flexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))==3,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[3])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")])
,0) )
flexible.detect <- mean(flexible.detect)
moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:7)]) &
s[x,"real_dif_2"]%in%c(s[x,paste0("dif_detect_",1:7)]) & s[x,"real_dif_3"]%in%c(s[x,paste0("dif_detect_",1:7)]),0) )
moreflexible.detect <- mean(moreflexible.detect)
}
if (nb.dif.true==1 & unique(b$J)==4 & unique(b$M)==2) {
perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))==1, s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1")])
,0) )
prop.perfect <- mean(perfect.detection)
flexible.detect <- prop.perfect
moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:4)]),0) )
moreflexible.detect <- mean(moreflexible.detect)
percent.detect <- sapply(1:1000,function(x) sum(s[x,paste0("dif_detect_",1:4)]%in%c(s[x,c("real_dif_1")]))/1)
percent.detect <- mean(percent.detect)
}
if (nb.dif.true==2 & unique(b$J)==4 & unique(b$M)==2) {
perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))==2,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:4),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2")])
,0) )
prop.perfect <- mean(perfect.detection)
flexible.detect <- prop.perfect
moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:4)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:4)]) &
s[x,"real_dif_2"]%in%c(s[x,paste0("dif_detect_",1:4)]),0) )
moreflexible.detect <- mean(moreflexible.detect)
percent.detect <- sapply(1:1000,function(x) sum(s[x,paste0("dif_detect_",1:4)]%in%c(s[x,c("real_dif_1","real_dif_2")]))/2)
percent.detect <- mean(percent.detect)
}
if (nb.dif.true==2 & unique(b$J)==7 & unique(b$M)==2) {
perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))==2,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2")])
,0) )
prop.perfect <- mean(perfect.detection)
flexible.detect <- prop.perfect
moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:7)]) &
s[x,"real_dif_2"]%in%c(s[x,paste0("dif_detect_",1:7)]),0) )
moreflexible.detect <- mean(moreflexible.detect)
percent.detect <- sapply(1:1000,function(x) sum(s[x,paste0("dif_detect_",1:7)]%in%c(s[x,c("real_dif_1","real_dif_2")]))/2)
percent.detect <- mean(percent.detect)
}
if (nb.dif.true==3 & unique(b$J)==7 & unique(b$M)==2) {
perfect.detection <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))==3,s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[1])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[2])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]) & s[x,paste0('dif_detect_',which(sapply(paste0("dif_detect_",1:7),function(y) !is.na(s[x,y])))[3])]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")])
,0) )
prop.perfect <- mean(perfect.detection)
flexible.detect <- prop.perfect
moreflexible.detect <- sapply(1:1000,function(x) ifelse(sum(!is.na(s[x,paste0("dif_detect_",1:7)]))!=0,s[x,"real_dif_1"]%in%c(s[x,paste0("dif_detect_",1:7)]) &
s[x,"real_dif_2"]%in%c(s[x,paste0("dif_detect_",1:7)]) & s[x,"real_dif_3"]%in%c(s[x,paste0("dif_detect_",1:7)]),0) )
moreflexible.detect <- mean(moreflexible.detect)
percent.detect <- sapply(1:1000,function(x) sum(s[x,paste0("dif_detect_",1:7)]%in%c(s[x,c("real_dif_1","real_dif_2","real_dif_3")]))/3)
percent.detect <- mean(percent.detect)
}
z <- data.frame(m.beta=mean(s$beta),
se.empirical.beta=sd(s$beta),
se.analytical.beta=mean(s$se_beta),
@ -693,7 +802,12 @@ compile_simulation2_resali <- function(scenario) {
true.value.in.ci.p=mean(true.value.in.ci),
h0.rejected.p=mean( (s$beta-1.96*s$se_beta)>0 | (s$beta+1.96*s$se_beta)<0 ),
beta.same.sign.truebeta.p=mean(beta.same.sign.truebeta.p),
beta.same.sign.truebeta.signif.p=mean(beta.same.sign.truebeta.p[num.reject])
beta.same.sign.truebeta.signif.p=mean(beta.same.sign.truebeta.p[num.reject]),
dif.detected=dif.d,
prop.perfect=prop.perfect,
flexible.detect=flexible.detect,
moreflexible.detect=moreflexible.detect,
percent.detect=ifelse(name%%2==0,NA,percent.detect)
)
d <- cbind(b,a,z)
d$prop.

8
RProject/Scripts/functions/resali.R
Unescape Escape View File

@ -48,8 +48,8 @@ resali <- function(df=NULL,items=NULL,group=NULL,verbose=T) {
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)]
dat[dat$TT==1,paste0("item",res.item,'TT')] <- dat[dat$TT==1,paste0('item',res.item)]
dat[dat$TT==0,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)
@ -118,8 +118,8 @@ resali <- function(df=NULL,items=NULL,group=NULL,verbose=T) {
res.item <- gsub("[a-z]", "",colnames(resp)[which.max(fval)])
res.items <- c(res.items,res.item)
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)]
dat[dat$TT==1,paste0("item",res.item,'TT')] <- dat[dat$TT==1,paste0('item',res.item)]
dat[dat$TT==0,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)

Write Preview
Loading…
Cancel
Save