Rasch-dikotomis

# install.packages("eRm")
# install.packages("ltm")
# install.packages("difR")
# install.packages("MASS")
# install.packages("xtable")

library("eRm")
library("ltm")
library("difR")
library("xtable")

#Pastikan nama file CSV sesuai dengan nama file yang dimiliki
#Kebetulan data kita kolom paling atas bernama Person, nantinya desesiaukan dengan penamaan yang kita berikan
raw <- read.csv(file = 'NilaiTes.csv', row.names = "Person")
sink('log/output.txt')
# Data skor adalah sepuluh kolom pertama setelah person, maka 1:10
score <- raw[,1:24]
print("Data Score: ")
print(score)


# Menyimpan data log, pastikan nama folder disesuaikan dengan
# memngekspor data ke latex, kalau suatu saat dibutuhkan
print(xtable(score, type = "latex"), file = "log/latex.tex")

# Menginisiasi Rasch model untuk 1 parameter logistik (1PL)
res_rm_1 <- RM(score)
# Mencetak data
print("Hasil iterasi dengan Rasg Model 1Pl : ")
print(res_rm_1)
# Mencetak ringkasan data
print("Ringkasan data : ")
summary(res_rm_1)

# Menghitung nilai indek kesukaran (difficulty logit)
print("Indek kesukaran butir : ")
betas <- -coef(res_rm_1)

# mencetak grafik difficulty logit item no 1-12
# Jika aitemnya bnayk, untuk memudahkan analisis bisa dicetak per-12 item
png(file="plot/ICC_1_12.png", width=600, height=350)
# 1:10 artinya item no 1 - 12
plotjointICC(res_rm_1, item.subset =  1:12, cex = .6)
# Menambahkan sumbu vertikal berwarna abu-abu, silahkan ganti warna sesuai selera
abline(v = -0.18, col = "grey")
# Menambahkan sumbu horizontal berwarna abu-abu, silahkan ganti warna sesuai selera
abline(h = .5, col = "grey")
# Menambahkan sumbu memabtu kita untuk memperlihatkan item fit dan difficulty logit dari setiap item,
# Jika item bersifat normal, berupa garis menyerupai huruf S, dari kiri ke kanan naik
dev.off()

# mencetak grafik difficulty logit item no 13-24
png(file="plot/ICC_13_24.png", width=600, height=350)
# 1:10 artinya item no 13 - 24
plotjointICC(res_rm_1, item.subset =  13:24, cex = .6)
# Menambahkan sumbu vertikal berwarna abu-abu, silahkan ganti warna sesuai selera
abline(v = -0.18, col = "grey")
# Menambahkan sumbu horizontal berwarna abu-abu, silahkan ganti warna sesuai selera
abline(h = .5, col = "grey")
# Menambahkan sumbu memabtu kita untuk memperlihatkan item fit dan difficulty logit dari setiap item,
# Jika item bersifat normal, berupa garis menyerupai huruf S, dari kiri ke kanan naik
dev.off()

# Mengurutkan nilai indek kesukaran
print("Ringkasan Indek kesukaran : ")
round(sort(betas), 5)
# Menampilkan rata-rata dari indeks kesukaran
print("Nilai rerata: ")
mean(betas[ 1:24])           # Want-items

# Menampilkan peta  person-item, dari plot ini kita bisa mengetahui item maupun person yang non fit
# Item non fit adalah item yang tidak mengukur sesuai dengan indikator
# Person non fit adalah person yang tidak sungguh-sungguh dalam pengukuran
png(file="plot/Map Person-item.png",width=600, height=350)
plotPImap(res_rm_1, cex.gen = .55)
dev.off()

# Map person - item dirutkan, untuk memudahkan
png(file="plot/Map Person-item-sorted.png",width=600, height=350)
plotPImap(res_rm_1, cex.gen = .55, sorted = TRUE)
dev.off()

# Menampilkan data normal
tmp1 <- RM(score, sum0 = FALSE)
print("Menampilkan data ability  : ")
print(tmp1)
print("Rerata ability person  : ")
round(coef(tmp1), 5)


# Menginisiasi Rasch model untuk 2 parameter logistik (2PL)
res_rm_2 <- rasch(score)
print("Analisis Rasch model 2PL  : ")
print(res_rm_2)

# Menghitung korelasi antara koefisien logistik 2 PL dengan 1 PL
print("Korelasi antara 2PL dengan 1PL  : ")
cor(coef(res_rm_2)[, 1], betas)
# A 2PL model cannot be estimated with eRm, but with the package ltm, for example.
# ltm() takes a formula with a tilde (~) as its first argument. The left-hand side
# must be the data and the right-hand side are the latent variables (only one
# here).


res_2pl_1 <- ltm(score ~ z1)
print("Menentukan nilai estimasi untuk 2PL : ")
print(res_2pl_1)

# Plot karakteristik item
png(file="plot/Item_characteristic_1_24.png",width=600, height=350)
plot(res_2pl_1, items =  1:24)
abline(v = -0.18, col = "grey")
abline(h = .5, col = "grey")
dev.off()

# Menentukan nilai analisis ANOVA antara 1 PL dengan 2 PL

print("Hasil analsisi anova 1PL dengan 2PL : ")
anova(res_rm_2, res_2pl_1)
# Korelasi antara 1PL dengan 2PL
print("Koefisien korelasi antara 1PL dengan 2PL : ")
cor(coef(res_rm_2)[, 1], coef(res_2pl_1)[, 1])


# Mereduksi data yang tidak fit

lrt_1 <- LRtest(res_rm_1, splitcr = raw$Gender)
png(file="plot/Graphical_model_check.png",width=600, height=350)
plotGOF(lrt_1, conf = list(), tlab = "number", xlab = "Women", ylab = "Men")
dev.off()
print(lrt_1)

tmp1 <- RM(score[, -6])
print("Data setelah direduksi: ")
LRtest(tmp1, splitcr = raw$Gender)
print(tmp1)

# Uji Wald
Waldtest(res_rm_1, splitcr = raw$Gender)

# Uji Item fit
pp_ml_1 <- person.parameter(res_rm_1)
itemfit(pp_ml_1)

# Plot item map
png(file="plot/ItemMap.png",width=600, height=350)
plotPWmap(res_rm_1)
dev.off()

# Uji Mantel-Haenszel
# Detection of Differential Item Functioning
tmp1 <- difMH(score, group = raw$Gender, focal.name = 1)
tmp1

png(file="plot/Uji mantel-Hanzel.png", width=600, height=350)
plot(tmp1)
dev.off()

# Pengujian cgi-square untuk uji Lord
print("Nilai Chi-Square Uji Lord 1 PL")
tmp1 <- difLord(score, group = raw$Gender, focal.name = 1, model = "1PL", discr = NULL)

# Kurva hasil uji chi-square
png(file="plot/Chi-square Pengaruh Gender.png", width=600, height=350)
plot(tmp1)
dev.off()

# Data plot individual, contohnya no. 6
png(file="plot/Item_no_6.png",
width=600, height=350)
plot(tmp1, plot = "itemCurve", item = 6)
dev.off()

# Parameter person
# Nilai esrimasi maksimum-likehood untuk person (ML)

print("Nilai estimasi maksimum-likehood 1PL: ")
tmp1 <- person.parameter(res_rm_1)
pp_ml <- coef(tmp1)
# MAP dan EAP
pp_map <- factor.scores(res_rm_2, method = "EB", resp.patterns = score)
pp_eap <- factor.scores(res_rm_2, method = "EAP", resp.patterns = score)
tmp1 <- data.frame(ML = pp_ml, MAP = pp_map$score.dat$z1, EAP = pp_eap$score.dat$z1)
round(cor(tmp1), 4)

png(file="plot/Kurva_ML.png", width=600, height=350)
plot(tmp1[, 1:2])
dev.off()


print("Nilai estimasi maksimum-likehood 2PL: ")
pp_2pl <- factor.scores(res_2pl_1, method = "EB", resp.patterns = score)
cor(pp_map$score.dat$z1, pp_2pl$score.dat$z1)
#> [1] 0.9956129

# Nilai item tes no 1-10
res_2pl_1
png(file="plot/Item_test_1_10.png",
width=600, height=350)
plot(res_2pl_1, items =  1:10, type = "IIC", ylim = c(0, 1.3))
dev.off()

# Kurva standard
png(file="plot/Item_test_standar.png",
width=600, height=350)
plot(res_2pl_1, items = 0, type = "IIC")
dev.off()
sink()