library(dplyr)
library(tidyr)
library(lavaan)
library(mirt)
library(psych)
library(lavaanPlot)
library(ggplot2)

# Binary Data

Adjustment <- read.csv("https://bit.ly/4kqbOcd")
head(Adjustment)

# Rename items for clarity
Adjustment <- Adjustment %>%
  dplyr::rename(Failure = I01, FallApart = I02, Wondering = I03, 
                HoldGrudges = I04, Hopeless = I05)
head(Adjustment)

describe(Adjustment[ , 2:6])
itemstats(Adjustment[ , 2:6])

tetrachoric(Adjustment[ , 2:6]) 
cor(Adjustment[ , 2:6]) # Pearson correlations are attenuated

AdjustmentItems <- select(Adjustment, c("Failure", "FallApart", "Wondering", "HoldGrudges", "Hopeless"))

# Specify the model
Adjustment.OneF <- 'Adjustment =~ Failure + FallApart + Wondering + HoldGrudges + Hopeless'

# Fit the model
Adjustment.OneF.Fit <- cfa(Adjustment.OneF, data = AdjustmentItems,
                           ordered = c("Failure", "FallApart", "Wondering", "HoldGrudges", "Hopeless"),
                           std.lv = TRUE)

summary(Adjustment.OneF.Fit, fit.measures = TRUE)

lavaanPlot(model = Adjustment.OneF.Fit,
           node_options = list(shape = "box", fontname = "Helvetica"),
           edge_options = list(color = "grey"), coefs = T)

TwoPL <- mirt(AdjustmentItems, model = 1, SE = TRUE, verbose = F)
TwoPL

coef(TwoPL, IRTpars = TRUE, simplify = TRUE)
coef(TwoPL, IRTpars = TRUE, printSE = TRUE)
coef(TwoPL, simplify = TRUE)

plot(TwoPL, type = "trace")
plot(TwoPL, type = "trace", facet = FALSE)

OnePLModel <- mirt.model('Theta = 1-5
                    CONSTRAIN = (1-5, a1)')
OnePL <- mirt(AdjustmentItems, OnePLModel, verbose = F)
OnePL

coef(OnePL, IRTpars = TRUE, simplify = TRUE)
plot(OnePL, type = "trace", facet = FALSE)
anova(OnePL, TwoPL)

M2(TwoPL)
itemfit(TwoPL, fit_stats = "S_X2", p.adjust = "BH")
residuals(TwoPL, type = "Q3")
residuals(TwoPL, type = "LD")

# Stanardize the X2-LD statistics
LD <- residuals(TwoPL, type = "LD")
k <- 2 # Number of categories for one item
m <- 2 # Number of categories for the other item
(abs(LD)-(k-1)*(m-1))/sqrt(2*(k-1)*(m-1))

plot(TwoPL, type = 'infotrace', facet_items = FALSE)

plot(TwoPL, type = "info")
plot(TwoPL, type = "infoSE")
plot(TwoPL, type = "itemscore")
plot(TwoPL, type = "score")

fscores(TwoPL, method = "EAP", full.scores = FALSE)

fscores(TwoPL, method = "MAP", full.scores = FALSE)

AdjustmentEAP <- fscores(TwoPL, method = "EAP", full.scores = TRUE, full.scores.SE = TRUE) %>%
  as.data.frame() %>%
  dplyr::rename(Theta_EAP = F1, SE_EAP = SE_F1)  # Rename columns for clarity 

Adjustment <- Adjustment %>%
  mutate(Theta_EAP = AdjustmentEAP$Theta_EAP,
         SE_EAP = AdjustmentEAP$SE_EAP)

head(Adjustment)



# Polytomous Data

MFW <- read.csv("https://bit.ly/4hfjZFo")

MFW <- MFW %>%
  dplyr::rename(ColdWeather = mfw1, VitaminC = mfw2, ChickenSoup = mfw3,
                WashingHands = mfw6, MultivitaminsCold = mfw7, CarbonatedBev = mfw8,
                PregnantPeriod = mfw9, WhiteSpots = mfw10, ShowerAfterSex = mfw11,
                Knuckles = mfw12, StarveFever = mfw13)
head(MFW)

GRM <- mirt(MFW, model = 1,
            SE = T, verbose = F)
GRM

coef(GRM, IRTpars = T, simplify = T)

M2(GRM, type = "C2")
itemfit(GRM, na.rm = TRUE, p.adjust = "BH")

LD <- residuals(GRM)

k <- 4 # Number of categories for one item
m <- 4 # Number of categories for the other item
(abs(LD)-(k-1)*(m-1))/sqrt(2*(k-1)*(m-1))

plot(GRM, type = "trace")
plot(GRM, type = "infotrace")

plot(GRM, type = "infoSE")

MFW_Theta <- fscores(GRM, method = "EAP", full.scores = TRUE, full.scores.SE = TRUE)
MFW_Theta <- data.frame(MFW_Theta)
colnames(MFW_Theta) <- c("GRM_EAP", "GRM_SE")
MFW <- cbind(MFW, MFW_Theta)
head(MFW)

ggplot(MFW, aes(x = GRM_EAP)) +
  geom_histogram(bins = 30, fill = "skyblue", color = "black") +
  xlab("IRT EAPs (Medical Folk Wisdom)") +
  ylab("Frequency") +
  theme_minimal(base_size = 14)

MFW$Sum_Score <- rowSums(MFW, na.rm = TRUE)
ggplot(MFW, aes(x = Sum_Score, y = GRM_EAP)) +
  geom_point() +
  theme_minimal(base_size = 12)

GPCM <- mirt(MFW[1:11], model = 1,
             itemtype = "gpcm", verbose = F)
coef(GPCM, IRTpars = T, simplify = T)

plot(GPCM, type = "trace", which.item = 9)

GPCM <- mirt(MFW[1:11], model = 1,
             itemtype = "gpcm", verbose = F)
coef(GPCM, IRTpars = T, simplify = T)

plot(GPCM, type = "trace", which.item = 4)
plot(GPCM, type = "trace", facet = T)

NRM <- mirt(MFW[1:11], itemtype = "nominal", verbose = F)
coef(NRM, IRTpars = T, simplify = T)
NRM

anova(GRM, NRM)

plot(NRM, type = "trace")

NRMprob <- function(x) {
  exp(x) / sum(exp(x))
}

# Compute response probabilities for an item at theta = 0
theta <- 0  # Latent trait level
a_params <- c(-1.245, -0.573, 0.057, 1.761)  # Discrimination parameters
c_params <- c(1.981, 0.724, -0.286, -2.420)  # Intercept parameters

# Compute exponentiated values
exponents <- a_params * (theta - c_params)

# Compute and print probabilities
probabilities <- NRMprob(exponents)
probabilities

RWA <- read.csv("https://bit.ly/4ia6nws")
names(RWA)
RWA <- RWA %>%
  dplyr::rename(DefyAuthority = RWA_1, Discipline = RWA_2, Obedience = RWA_4,
                RespectAuthority = RWA_5, Protest = RWA_6, StraightNarrow = RWA_7,
                OldFashioned = RWA_8, GodsLaws = RWA_9, NudistCamps = RWA_10,
                YoungPeopleExp = RWA_11, PreMaritalSex = RWA_13, StrongGovt = RWA_14,
                KindCriminals = RWA_15, NoTougherGovt = RWA_16,
                PrisonsDisgrace = RWA_18, StrongMedicine = RWA_19,
                Age = agebinary)

RWA.GRM <- mirt(RWA[1:16], 1, itemtype = "graded", verbose = F)
coef(RWA.GRM, IRTpars = T, simplify = T)
M2(RWA.GRM, type = "C2")

plot(RWA.GRM, type = "trace", facet = TRUE)

RWA.binary <- RWA %>%
  dplyr::mutate(across(c(-Age), ~ case_when(
    . == 1 ~ 1,
    . == 2 ~ 2,
    . == 3 ~ 2,
    TRUE ~ . 
  )))

RWA.TwoPL <- mirt(RWA.binary[c(1:16)], 1, itemtype = "2PL", verbose = F)
coef(RWA.TwoPL, IRTpars = T, simplify = T)
M2(RWA.TwoPL)

plot(RWA.TwoPL, type = "trace", facet = TRUE)

RWA.binary$Age <- factor(RWA.binary$Age, labels = c("40 & Under", "Over 40"))
MultiGroupMod <- multipleGroup(RWA.binary[, c(1:16)], 1, group = RWA.binary$Age, SE = TRUE, verbose = F)
coef(MultiGroupMod, IRTpars = T, simplify = T)

a_DIF <- DIF(MultiGroupMod, 'a1', p.adjust = "BH")
a_DIF

b_DIF <- DIF(MultiGroupMod, c('d'), p.adjust = 'BH')
b_DIF

plot(MultiGroupMod, type = "trace", which.item = c(1, 7, 8, 10, 11))
plot(MultiGroupMod, which.item = c(1, 7, 8, 10, 11))

IgnoreDIF <- multipleGroup(RWA.binary[, c(1:16)], 1,
                           group = RWA.binary$Age,
                           invariance = c("slopes", "intercepts", "free_mean", "free_variance"),
                           SE = TRUE, verbose = F)
coef(IgnoreDIF, IRTpars = T, simplify = T)

par_table <- mod2values(MultiGroupMod)
par_table[ , c(1, 2, 4, 5)]

DIFfree <- list(c(6, 72), c(10, 76), c(14, 80), c(18, 84), 
                c(22, 88), c(34, 100), c(46, 112), c(50, 116), 
                c(54, 120), c(58, 124), c(62, 128))

AccountForDIF <- multipleGroup(RWA.binary[, c(1:16)], 1,
                               group = RWA.binary$Age, constrain = DIFfree,
                               invariance = c("slopes", "free_mean", "free_var"),
                               SE = TRUE, verbose = F)

coef(AccountForDIF, IRTpars = T, simplify = T)