####################################################################################
#                        L9 - DIFFERENTIAL ITEM FUNCTIONING                        #
#                                    NMST570                                       #
#                             Last update: 10/12/2019                              #
####################################################################################

# Slides available at http://www.cs.cas.cz/martinkova/NMST570.html

####################################################################################
# Corrections for multiple comparison - toy example ################################
####################################################################################

# ordered p-value
pval <- c(0.001, 0.004, 0.006, 0.018, 0.021, 0.031, 0.042, 0.243, 0.362, 0.783)
# significance level
alpha <- 0.05
# number of items
I <- length(pval)

# Bonferroni correction
### boundary
alpha/I
### index k
(k_bonf <- which(pval > alpha/I)[1])
### DIF items
c(1:10)[1:(k_bonf - 1)]
### adjusted p-value, compare with alpha
pmin(1, pval*I)
### with p.adjust function
p.adjust(pval, method = "bonferroni")

# Holm correction
### boundary
alpha/(I + 1 - 1:I)
### index k
(k_holm <- which(pval > alpha/(I + 1 - 1:I))[1])
### DIF items
c(1:10)[1:(k_holm - 1)]
### adjusted p-value, compare with alpha
pmin(1, cummax(pval*(I + 1 - 1:I)))
### with p.adjust function
p.adjust(pval, method = "holm")

# Benjamini-Hochberg correction
### boundary
1:I * alpha/I
### index k
(k_holm <- rev(which(pval <= 1:I * alpha/I))[1])
### DIF items
c(1:10)[1:(k_holm)]
### adjusted p-value, compare with alpha
rev(pmin(1, cummin(I/I:1 * rev(pval))))
### with p.adjust function
p.adjust(pval, method = "BH")