# Project: The 19-Item Environmental Knowledge Test (EKT-19): A Short, Psychometrically Robust Measure of Environmental Knowledge
# Last Edit: 27/09/22
# Analysis: IRT

# Install IRT packages
install.packages("ltm")
install.packages("dplyr")
install.packages("readxl")

# Read in libraries
library(ltm)
library(dplyr)
library(readxl)

# Begin saving data output in separate file
sink(file = "EK_IRT_Output.txt") 

# Import data
EK_irt <- read_excel("C:\\Users\\lkp32\\R\\EnvKnowledge\\EKT-19_Cleaned.xlsx")
View(EK_irt) # Visually inspect data

# Remove failed attention checks
EK_irt_filtered <- filter(EK_irt, Attention == 1)
View(EK_irt_filtered)

# Filter to only include EK items in analysis
EK_irt_EKItemsOnly <- select(EK_irt_filtered, EK1:EK15, EK16:EK22, EK23:EK30)
class(1:30) # Check variable type
View(EK_irt_EKItemsOnly)

# Creating a Rasch (1PL) Model 
IRTrasch <- rasch(EK_irt_EKItemsOnly, IRT.param = TRUE, constraint = cbind(length(EK_irt_EKItemsOnly) + 1, 1))
print(IRTrasch)
summary(IRTrasch)
coef(IRTrasch)

# Creating a dichotomous IRT (2PL)
IRTmodel <- ltm(EK_irt_EKItemsOnly ~ z1, IRT.param = TRUE)
print(IRTmodel)
summary(IRTmodel)
coef(IRTmodel)

# Testing assumption of unidimensionality using Modified Parallel Analysis
unidim <- unidimTest(IRTmodel, EK_irt_EKItemsOnly)
unidim
plot(unidim, type = "b", pch = 1:2)
legend("topright", c("Real Data", "Average Simulated Data"), lty = 1, 
       pch = 1:2, col = 1:2, bty = "n")

# Calculating ability estimates for each person
est <- factor.scores(IRTmodel, method = "MI", B= 20)
options(max.print = 100000)
est
est$score.dat$z1[1:5]
mean(est$score.dat$z1)

# Getting parameters for each model
coef(IRTrasch, prob = TRUE)
coef(IRTmodel, prob = TRUE)

# Comparing 1PL and 2PL models
anova(IRTrasch, IRTmodel)

# Plotting IRTmodel
plot(IRTmodel, type = "ICC") # Plot all items together
plot(IRTmodel, type = "ICC", items = 19) # Plot one item at a time, by item number
plot(IRTmodel, type = "ICC", items = c(2,3)) # Plot several items at a time, by item number

# Set colour palette for plots
palette("Set2")

# Creating item information curves for each item (3 per graph) 
plot(IRTmodel, type = "IIC", items = 1:3, lwd = 2, cex.lab = 1.2, col = palette(), xlab = "Environmental Knowledge", ylab = "Information", main = NULL, labels = NULL, legend = TRUE, cx = "topright", cy = NULL, ncol = 1, bty = "n") 
plot(IRTmodel, type = "IIC", items = 4:6, lwd = 2, cex.lab = 1.2, col = palette(), xlab = "Environmental Knowledge", ylab = "Information", main = NULL, labels = NULL, legend = TRUE, cx = "topright", cy = NULL, ncol = 1, bty = "n")
plot(IRTmodel, type = "IIC", items = 7:9, lwd = 2, cex.lab = 1.2, col = palette(), xlab = "Environmental Knowledge", ylab = "Information", main = NULL, labels = NULL, legend = TRUE, cx = "topright", cy = NULL, ncol = 1, bty = "n") 
plot(IRTmodel, type = "IIC", items = 10:12, lwd = 2, cex.lab = 1.2, col = palette(), xlab = "Environmental Knowledge", ylab = "Information", main = NULL, labels = NULL, legend = TRUE, cx = "topright", cy = NULL, ncol = 1, bty = "n")
plot(IRTmodel, type = "IIC", items = 13:15, lwd = 2, cex.lab = 1.2, col = palette(), xlab = "Environmental Knowledge", ylab = "Information", main = NULL, labels = NULL, legend = TRUE, cx = "topright", cy = NULL, ncol = 1, bty = "n") 
plot(IRTmodel, type = "IIC", items = 16:18, lwd = 2, cex.lab = 1.2, col = palette(), xlab = "Environmental Knowledge", ylab = "Information", main = NULL, labels = NULL, legend = TRUE, cx = "topright", cy = NULL, ncol = 1, bty = "n") 
plot(IRTmodel, type = "IIC", items = 19:21, lwd = 2, cex.lab = 1.2, col = palette(), xlab = "Environmental Knowledge", ylab = "Information", main = NULL, labels = NULL, legend = TRUE, cx =  "topleft", cy = NULL, ncol = 1, bty = "n") 
plot(IRTmodel, type = "IIC", items = 22:24, lwd = 2, cex.lab = 1.2, col = palette(), xlab = "Environmental Knowledge", ylab = "Information",main = NULL, labels = NULL, legend = TRUE, cx = "topright", cy = NULL, ncol = 1, bty = "n") 
plot(IRTmodel, type = "IIC", items = 25:27, lwd = 2, cex.lab = 1.2, col = palette(), xlab = "Environmental Knowledge", ylab = "Information", main = NULL, labels = NULL, legend = TRUE, cx = "topright", cy = NULL, ncol = 1, bty = "n") 
plot(IRTmodel, type = "IIC", items = 28:30, lwd = 2, cex.lab = 1.2, col = palette(), xlab = "Environmental Knowledge", ylab = "Information", main = NULL, labels = NULL, legend = TRUE, cx = "topright", cy = NULL, ncol = 1, bty = "n") 

# Generates pattern of data, expected scores for that person
person.fit(IRTmodel) 
item.fit(IRTmodel)

# Creating 3PL model, which includes guessing dimension 
IRTmodel2 = tpm(EK_irt_EKItemsOnly, type = "latent.trait", IRT.param = TRUE)
summary(IRTmodel2)
coef(IRTmodel2)
plot(IRTmodel2, type = "ICC") # Plot IIC curve
factor.scores(IRTmodel2)
person.fit(IRTmodel2)
item.fit(IRTmodel2)

# Comparing models to see if adding a guessing paramater improved model
anova(IRTmodel, IRTmodel2)

# saving output file
sink(file = NULL)