resp <- df[,c(dif.items.list)]
}
}
else {
resp <- df[,c(nodif.items.list)]
}
print(grp)
tam1 <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
}
}
return(summary(tam1))
}
rosali(dat=dat,items=1:7,group="TT")
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
return(c(model_a,model_b))
}
rosali(dat=dat,items=1:7,group="TT")
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])-1
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
return(c(model_a,model_b))
}
rosali(dat=dat,items=1:7,group="TT")
## File Name: pcm.R
## File version: 1.0
#' @import TAM
#' @export
pcm <- function(df=NULL,items=NULL,group=NULL,model="PCM2",method="MML",dif.items=NULL) {
# Prepare analysis
if (is.null(items)) {
nbitems <- sum(sapply(1:100,function(x) paste0('item',x)) %in% colnames(df))
items <- paste0('item',seq(1,nbitems))
resp <- df[,items]
}
else {
nbitems <- length(items)
resp <- df[,paste0("item",items)]
}
if (is.null(group)) {
grp <- NULL
}
else {
grp <- df[,group]
df$grp <- grp
if (!is.null(dif.items)) {
for (i in dif.items) {
df[,paste0('item',i,'_noTT')] <- NA
df[,paste0('item',i,'_TT')] <- NA
df[df$grp==0,paste0('item',i,'_noTT')] <- df[df$grp==0,paste0('item',i)]
df[df$grp==1,paste0('item',i,'_TT')] <- df[df$grp==1,paste0('item',i)]
}
}
}
# Analyze
if (method=='MML') {
if (is.null(dif.items)) {
tam1 <-  TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
}
else {
nodif.items <- setdiff(items,dif.items)
nodif.items.list <- paste0("item",nodif.items)
dif.items.list <- c(paste0("item",dif.items,"_noTT"),paste0("item",dif.items,"_TT"))
if (!is.null(dif.items)) {
if (length(nodif.items)!=0) {
resp <- df[,c(nodif.items.list,dif.items.list)]
}
else {
resp <- df[,c(dif.items.list)]
}
}
else {
resp <- df[,c(nodif.items.list)]
}
tam1 <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
}
}
return(summary(tam1))
}
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])-1
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
return(c(model_a))
}
rosali(dat=dat,items=1:7,group="TT")
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])-1
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
return(c(model_b))
}
rosali(dat=dat,items=1:7,group="TT")
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])-1
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
}
rosali(dat=dat,items=1:7,group="TT")
model_b <- pcm(df=dat,items=1:7,group = "TT",dif.items = 1:7)
model_a <- pcm(df=dat,items=1:7,group = "TT",dif.items = NULL)
lrt
lrt''
??lrt
loglik(model_a)
logLik(model_a)
## File Name: pcm.R
## File version: 1.0
#' @import TAM
#' @export
pcm <- function(df=NULL,items=NULL,group=NULL,model="PCM2",method="MML",dif.items=NULL) {
# Prepare analysis
if (is.null(items)) {
nbitems <- sum(sapply(1:100,function(x) paste0('item',x)) %in% colnames(df))
items <- paste0('item',seq(1,nbitems))
resp <- df[,items]
}
else {
nbitems <- length(items)
resp <- df[,paste0("item",items)]
}
if (is.null(group)) {
grp <- NULL
}
else {
grp <- df[,group]
df$grp <- grp
if (!is.null(dif.items)) {
for (i in dif.items) {
df[,paste0('item',i,'_noTT')] <- NA
df[,paste0('item',i,'_TT')] <- NA
df[df$grp==0,paste0('item',i,'_noTT')] <- df[df$grp==0,paste0('item',i)]
df[df$grp==1,paste0('item',i,'_TT')] <- df[df$grp==1,paste0('item',i)]
}
}
}
# Analyze
if (method=='MML') {
if (is.null(dif.items)) {
tam1 <-  TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
}
else {
nodif.items <- setdiff(items,dif.items)
nodif.items.list <- paste0("item",nodif.items)
dif.items.list <- c(paste0("item",dif.items,"_noTT"),paste0("item",dif.items,"_TT"))
if (!is.null(dif.items)) {
if (length(nodif.items)!=0) {
resp <- df[,c(nodif.items.list,dif.items.list)]
}
else {
resp <- df[,c(dif.items.list)]
}
}
else {
resp <- df[,c(nodif.items.list)]
}
tam1 <- TAM::tam.mml(resp=resp,Y=grp,irtmodel = model,est.variance = T,verbose=F)
}
}
return(tam1)
}
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])-1
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
}
model_b <- pcm(df=dat,items=1:7,group = "TT",dif.items = 1:7)
model_a <- pcm(df=dat,items=1:7,group = "TT",dif.items = NULL)
logLik(model_a)
logLik(model_b)
library(lmtest)
lrtest(model_a,model_b)
lrtest(model_a,model_b)$p
ss <- lrtest(model_a,model_b)
ss$`Pr(>Chisq)`
ss$`Pr(>Chisq)`[2]
a <- c()
a
c(a,1)
which(1:4,3)
?which
which.min(1:4)
which.min(2:4)
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])-1
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
ltest <- lmtest::lrtest(model_a,model_b)
if (ltest$`Pr(>Chisq)`[2]<1) {
model_c <- model_b
difit <- c()
a <- sapply(1:nbitems,function(x)
lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
while (min(a)<0.4) {
difit <- c(difit,which.min(a))
model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
a <- sapply(1:nbitems,function(x)
lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
}
}
}
rosali(dat=dat,items=1:7,group="TT")
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])-1
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
ltest <- lmtest::lrtest(model_a,model_b)
if (ltest$`Pr(>Chisq)`[2]<1) {
model_c <- model_b
difit <- c()
a <- sapply(1:nbitems,function(x)
lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
while (min(a)<0.4) {
difit <- c(difit,which.min(a))
model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
print(difit)
a <- sapply(1:nbitems,function(x)
lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
}
}
}
rosali(dat=dat,items=1:7,group="TT")
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])-1
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
ltest <- lmtest::lrtest(model_a,model_b)
if (ltest$`Pr(>Chisq)`[2]<1) {
model_c <- model_b
difit <- c()
a <- sapply(1:nbitems,function(x)
lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
while (min(a)<0.4) {
difit <- c(difit,which.min(a))
model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
print(difit)
a <- sapply(1:nbitems,function(x)
lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
print(a)
}
}
}
rosali(dat=dat,items=1:7,group="TT")
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])-1
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
ltest <- lmtest::lrtest(model_a,model_b)
if (ltest$`Pr(>Chisq)`[2]<1) {
model_c <- model_b
difit <- c()
a <- sapply(1:nbitems,function(x)
lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
while (min(a)<0.4) {
difit <- c(difit,which.min(a))
model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
print(difit)
a <- sapply(1:nbitems,function(x)
ifelse(x%in%difit,lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
print(a)
}
}
}
rosali(dat=dat,items=1:7,group="TT")
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])-1
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
ltest <- lmtest::lrtest(model_a,model_b)
if (ltest$`Pr(>Chisq)`[2]<1) {
model_c <- model_b
difit <- c()
a <- sapply(1:nbitems,function(x)
lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
while (min(a)<0.4) {
difit <- c(difit,which.min(a))
model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
print(difit)
a <- sapply(1:nbitems,function(x)
ifelse(!(x%in%difit),lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
print(a)
}
}
}
rosali(dat=dat,items=1:7,group="TT")
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])-1
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
ltest <- lmtest::lrtest(model_a,model_b)
difit <- c()
if (ltest$`Pr(>Chisq)`[2]<0.05) {
model_c <- model_b
a <- sapply(1:nbitems,function(x)
lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
while (min(a)<0.05) {
difit <- c(difit,which.min(a))
model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
print(difit)
a <- sapply(1:nbitems,function(x)
ifelse(!(x%in%difit),lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
}
}
return(difit)
}
rosali(dat=dat,items=1:7,group="TT")
dat <- read.csv("/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N200/scenario_20A_200.csv")
dat <- dat[dat$replication==2,]
rosali(dat=dat,items=1:7,group="TT")
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])-1
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
ltest <- lmtest::lrtest(model_a,model_b)
difit <- c()
if (ltest$`Pr(>Chisq)`[2]<0.05) {
model_c <- model_b
a <- sapply(1:nbitems,function(x)
lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
while (min(a)<0.05) {
difit <- c(difit,which.min(a))
model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
a <- sapply(1:nbitems,function(x)
ifelse(!(x%in%difit),lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
}
}
return(difit)
}
rosali(dat=dat,items=1:7,group="TT")
dat <- read.csv("/home/corentin/Documents/These/Recherche/Simulations/Data/DIF/N200/scenario_20A_200.csv")
dat <- dat[dat$replication==3,]
rosali(dat=dat,items=1:7,group="TT")
rosali <- function(dat=NULL,items=NULL,group=NULL) {
nbitems <- length(items)
items2 <- items
# Un seul groupe
if (length(group)!=1) {
stop('Only one variable can be used for the group option')
}
# Recoder groupe en 0/1
dat[,group] <- as.factor(dat[,group])
if (!(all(names(levels(dat[,group]))==c("0","1")))) {
levels(dat[,group]) <- c("0","1")
}
dat[,group] <- as.numeric(dat[,group])-1
model_a <- pcm(df=dat,items=items,group = group,dif.items = items)
model_b <- pcm(df=dat,items=items,group = group,dif.items = NULL)
ltest <- lmtest::lrtest(model_a,model_b)
difit <- c()
if (ltest$`Pr(>Chisq)`[2]<0.05) {
model_c <- model_b
a <- sapply(1:nbitems,function(x)
lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2])
k <- 0
while (min(a)<0.05/(nbitems-k)) {
difit <- c(difit,which.min(a))
model_c <- pcm(df=dat,items=items,group = group,dif.items = difit)
a <- sapply(1:nbitems,function(x)
ifelse(!(x%in%difit),lmtest::lrtest(model_c,pcm(df=dat,items=items,group = group,dif.items = c(difit,x)))$`Pr(>Chisq)`[2],1))
}
}
return(difit)
}
rosali(dat=dat,items=1:7,group="TT")
pcm(dat,items=1:7,group="TT",dif.items=1:2)
summary(pcm(dat,items=1:7,group="TT",dif.items=1:2))