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

Added custom ROSALI algorithm modifications

Browse Source
This commit is contained in:
Corentin Choisy
2024-03-04 16:58:39 +01:00
parent 22e6ac2dcd
commit ee41611941
6 changed files with 4646 additions and 133 deletions
Download Patch File Download Diff File Expand all files Collapse all files

85
Modules/ghquadm.ado
View File

@ -1,85 +0,0 @@
program define ghquadm
* stolen from gllamm6 who stole it from rfprobit (Bill Sribney)
version 4.0
parse "`*'", parse(" ")
local n = `1'
if `n' + 2 > _N {
di in red /*
*/ "`n' + 2 observations needed to compute quadrature points"
exit 2001
}
tempname x w xx ww a b
local i 1
local m = int((`n' + 1)/2)
matrix x = J(1,`m',0)
matrix w = x
while `i' <= `m' {
if `i' == 1 {
scalar `xx' = sqrt(2*`n'+1)-1.85575*(2*`n'+1)^(-1/6)
}
else if `i' == 2 { scalar `xx' = `xx'-1.14*`n'^0.426/`xx' }
else if `i' == 3 { scalar `xx' = 1.86*`xx'-0.86*x[1,1] }
else if `i' == 4 { scalar `xx' = 1.91*`xx'-0.91*x[1,2] }
else {
local im2 = `i' -2
scalar `xx' = 2*`xx'-x[1,`im2']
}
hermite `n' `xx' `ww'
matrix x[1,`i'] = `xx'
matrix w[1,`i'] = `ww'
local i = `i' + 1
}
if mod(`n', 2) == 1 { matrix x[1,`m'] = 0}
/* start in tails */
matrix `b' = (1,1)
matrix w = w#`b'
matrix w = w[1,1..`n']
matrix `b' = (1,-1)
matrix x = x#`b'
matrix x = x[1,1..`n']
/* other alternative (left to right) */
/*
above: matrix x = J(1,`n',0)
while ( `i'<=`n'){
matrix x[1, `i'] = -x[1, `n'+1-`i']
matrix w[1, `i'] = w[1, `n'+1-`i']
local i = `i' + 1
}
*/
matrix `2' = x
matrix `3' = w
end
program define hermite /* integer n, scalar x, scalar w */
* stolen from gllamm6 who stole it from rfprobit (Bill Sribney)
version 4.0
local n "`1'"
local x "`2'"
local w "`3'"
local last = `n' + 2
tempname i p
matrix `p' = J(1,`last',0)
scalar `i' = 1
while `i' <= 10 {
matrix `p'[1,1]=0
matrix `p'[1,2] = _pi^(-0.25)
local k = 3
while `k'<=`last'{
matrix `p'[1,`k'] = `x'*sqrt(2/(`k'-2))*`p'[1,`k'-1] /*
*/ - sqrt((`k'-3)/(`k'-2))*`p'[1,`k'-2]
local k = `k' + 1
}
scalar `w' = sqrt(2*`n')*`p'[1,`last'-1]
scalar `x' = `x' - `p'[1,`last']/`w'
if abs(`p'[1,`last']/`w') < 3e-14 {
scalar `w' = 2/(`w'*`w')
exit
}
scalar `i' = `i' + 1
}
di in red "hermite did not converge"
exit 499
end

1150
Modules/rosali_custom/rosali_nobf.ado Normal file
View File

File diff suppressed because it is too large Load Diff

1149
Modules/rosali_custom/rosali_nolrt.ado Normal file
View File

File diff suppressed because it is too large Load Diff

1149
Modules/rosali_custom/rosali_nolrt_nobf.ado Normal file
View File

File diff suppressed because it is too large Load Diff

1150
Modules/rosali_custom/rosali_original.ado Normal file
View File

File diff suppressed because it is too large Load Diff

96
RProject/.Rhistory
View File

@ -1,51 +1,3 @@
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])
)
d <- cbind(b,a,z)
d$prop.
return(d)
}
#### Compiled results
res.dat.dif <- compile_simulation2('5A_100')
for (x in results[seq(2,length(results))]) {
y <- compile_simulation2(x)
res.dat.dif <- bind_rows(res.dat.dif,y)
}
#### Create data.frame
results <- c(sapply(1:4,function(x) paste0(x,c('A','B','C','D','E'))),sapply(5:9,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
results2 <- c(sapply(10:20,function(x) paste0(x,c('A','B','C','D','E','F','G'))))
results <- c(sapply(c(50,100,200,300),function(x) paste0(results,'_',x)))
results2 <- c(sapply(c(50,100,200,300),function(x) paste0(results2,'_',x)))
results <- sort(results)
results2 <- sort(results2)
results <- c(results,results2)
#### Compiler function
compile_simulation <- 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<=4) {
s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N50/scenario_',scenario,'.csv'))
}
if (substr(scenario,start=nchar(scenario)-1,stop=nchar(scenario))=="50" & name>4) {
s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N50/scenario_',scenario,'.csv'))
}
if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name<=4) {
s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N100/scenario_',scenario,'.csv'))
}
if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="100" & name>4) {
s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N100/scenario_',scenario,'_nodif.csv'))
}
if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name<=4) {
s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N200/scenario_',scenario,'.csv'))
}
if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="200" & name>4) {
s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N200/scenario_',scenario,'_nodif.csv'))
}
if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name<=4) {
s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N300/scenario_',scenario,'.csv'))
}
if (substr(scenario,start=nchar(scenario)-2,stop=nchar(scenario))=="300" & name>4) {
s <- read.csv(paste0('/home/corentin/Documents/These/Recherche/Simulations/Analysis/NoDIF/N300/scenario_',scenario,'_nodif.csv'))
}
if (unique(s$J)==4) { if (unique(s$J)==4) {
if (unique(s$M)==2) { if (unique(s$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),m.item2=mean(s$item2),m.item3=mean(s$item3),m.item4=mean(s$item4))
@ -510,3 +462,51 @@ res.dat[res.dat$scenario %in% paste0(4,'D') & res.dat$N==50,]$theoretical.power
res.dat[res.dat$scenario %in% paste0(4,'E') & res.dat$N==50,]$theoretical.power <- 0.4328 res.dat[res.dat$scenario %in% paste0(4,'E') & res.dat$N==50,]$theoretical.power <- 0.4328
### DIF scenarios ### DIF scenarios
res.dat.dif$theoretical.power <- res.dat[81:nrow(res.dat),]$theoretical.power res.dat.dif$theoretical.power <- res.dat[81:nrow(res.dat),]$theoretical.power
install.packages("lordif")
library(lordif)
library(TAM)
library(doMC)
library(parallel)
library(pbmcapply)
library(funprog)
library(dplyr)
library(readxl)
dat1 <- read.csv(file = '/home/corentin/Documents/These/Recherche/Simulations/Data/NoDIF/N50/scenario_1A_50.csv')
dat1 <- dat1[dat1$replication==1,]
lordif(dat1,group=dat1$TT)
dat1
lordif(dat1,group=dat1$TT,selection = dat[,c('item1','item2',"item3",'item4')],model = "GPCM")
lordif(dat1,group=dat1$TT,selection = dat1[,c('item1','item2',"item3",'item4')],model = "GPCM")
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM")
dat1 <- read.csv(file = '/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N50/scenario_5A_50.csv')
dat1 <- dat1[dat1$replication==1,]
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM")
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000)
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta")
dat1 <- read.csv(file = '/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_5A_300.csv')
dat1 <- dat1[dat1$replication==1,]
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta")
dat1 <- read.csv(file = '/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_5D_300.csv')
dat1 <- dat1[dat1$replication==1,]
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta")
dat1 <- read.csv(file = '/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_7D_300.csv')
dat1 <- dat1[dat1$replication==1,]
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta")
dat1 <- dat1[dat1$replication<=20,]
dat1 <- read.csv(file = '/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_7D_300.csv')
dat1 <- dat1[dat1$replication==1,]
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta")
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GRM",maxIter = 10000,criterion = "Beta")
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta")
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta",MonteCarlo = T)
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta",inc=0.01)
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta",minCell=0.01)
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta",minCell=1)
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta",minCell=10)
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta")
dat1 <- read.csv(file = '/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N300/scenario_7A_300.csv')
dat1 <- dat1[dat1$replication<=20,]
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta")
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta",beta.change = 0.01)
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta",beta.change = 0.001)
lordif(dat1[,c('item1','item2',"item3",'item4')],group=dat1$TT,model = "GPCM",maxIter = 10000,criterion = "Beta",beta.change = 0.2)