####################################################################################
#                    L2 - RELIABILITY AND MEASUREMENT ERROR                        #
#                                    NMST570                                       #
#                             Last update: 13/10/2020                              #
####################################################################################

# Slides available at http://www.cs.cas.cz/martinkova/docsNMST570/NMST570_L2slides.pdf
# HW assignment available at http://www.cs.cas.cz/hladka/documents/NMST570_HW2_2020.pdf

####################################################################################
# Ex 3.1.1 #########################################################################
####################################################################################

# HCI_test_retest is a dataset closer described in McFarland et al. (2017)
# (see http://dx.doi.org/10.1187/cbe.16-10-0305)
# It includes answers from 45 students (`ID`) on 20 item test administered twice in row.
# variable `test` indicates whether results were obtained from 1st or 2nd administration.
# Dataset also includes information on obtained total score (`total`)

# load data
# you can use getwd() function to know where is your working directory
# you can change your working directory with function setwd(), see ?setwd
load("HCI_test_retest.rda")

# explore data
head(HCI_test_retest)
summary(HCI_test_retest)

# divide dataset by `test`
HCI_test   <- HCI_test_retest[HCI_test_retest$test == "1", ]
HCI_retest <- HCI_test_retest[HCI_test_retest$test == "2", ]

# sort datasets by `ID`
HCI_test   <- HCI_test[order(HCI_test$ID), ]
HCI_retest <- HCI_retest[order(HCI_retest$ID), ]

# check whether all IDs match in both datasets
all(HCI_test$ID == HCI_retest$ID)

# Pearson correlation coefficient for total scores
cor.test(HCI_test$total, HCI_retest$total)

####################################################################################
# Ex 3.1.2 #########################################################################
####################################################################################

# load packages
library(ShinyItemAnalysis)
library(psychometric)
library(psych)

# only item data
test <- HCI_test[, 3:22]; head(test)
retest <- HCI_retest[, 3:22]; head(retest)

#--------------------------------------------------------------
# 1. First-last half split approach
#--------------------------------------------------------------

# TEST
### splitting data, first 10 and last 10 items
test_first <- test[, 1:10]
test_last <- test[, 11:20]
### calculate total scores for subsets
test_score_first <- rowSums(test_first)
test_score_last <- rowSums(test_last)
### calculate correlation for total scores
cor(test_score_first, test_score_last)
### apply Spearmann-Brown formula with m = 2
2 * cor(test_score_first, test_score_last) / (1 + cor(test_score_first, test_score_last))

# TODO: Do the same for retest!
# HINT: In the previous code, change the test to retest.

#--------------------------------------------------------------
# 2. Odd and even half split approach
#--------------------------------------------------------------

# TEST
### splitting data, even and odd items
test_odd <- test[, seq(2, 20, 2)]
test_even <- test[, seq(1, 20, 2)]
### calculate total scores for subsets
test_score_odd <- rowSums(test_odd)
test_score_even <- rowSums(test_even)
### calculate correlation for total scores
cor(test_score_odd, test_score_even)
### apply Spearmann-Brown formula with m = 2
2 * cor(test_score_odd, test_score_even) / (1 + cor(test_score_odd, test_score_even))

# TODO: Do the same for retest!

#--------------------------------------------------------------
# 3. Random split
#--------------------------------------------------------------

# TEST
### splitting data, randomly choose 10 items from 20
set.seed(123)
samp <- sample(1:20, 10)
###### consider these columns into 1st sample
test_random1 <- test[, samp]
###### rest of columns into 2nd sample
test_random2 <- test[, setdiff(1:20, samp)]
### calculate total scores for subsets
test_score_random1 <- rowSums(test_random1)
test_score_random2 <- rowSums(test_random2)
### calculate correlation for total scores
cor(test_score_random1, test_score_random2)
### apply Spearmann-Brown formula with m = 2
2 * cor(test_score_random1, test_score_random2) / (1 + cor(test_score_random1, test_score_random2))

# TODO: Do the same for retest!

#--------------------------------------------------------------
# 4. Average of 10,000 split-halves
#--------------------------------------------------------------

# TEST
### splitting data, 10,000 random split-halves
set.seed(123)
test_split_10000 <- psych::splitHalf(test, raw = TRUE)
### calculate mean reliability
mean(test_split_10000$raw)

# TODO: Do the same for retest!

#--------------------------------------------------------------
# 5. Average of all possible split-halves
#--------------------------------------------------------------

# TEST
### splitting data, all possible split-halves
test_split_all <- psych::splitHalf(test, raw = TRUE, brute = TRUE)
### calculate mean reliability
mean(test_split_all$raw)

# TODO: Do the same for retest!

#--------------------------------------------------------------
# 6. The worst split
#--------------------------------------------------------------

# TEST
### calculate minimal reliability
min(test_split_all$raw)
# the same result
test_split_all$minrb

# TODO: Do the same for retest!

####################################################################################
# Ex 3.1.3 #########################################################################
####################################################################################

# TEST
### calculate Cronbach's alpha
(a_test <- psychometric::alpha(test))
### calculate CI
psychometric::alpha.CI(a_test, N = nrow(test), k = ncol(test), level = 0.95)

# TODO: Do the same for retest!

# You can also calculate Cronbach's alpha by hand:
### calculation of test score
test_score <- rowSums(test)
### variance of test score
(sigma2X <- var(test_score))
### variance of its items
(sigma2Xi <- sapply(test, var))
### number of items in test
(I <- ncol(test))
### formula for Cronbach's alpha:
(alpha_test <- I/(I - 1) * (1 - sum(sigma2Xi)/sigma2X))
# Compare with code of functions:
psychometric::alpha
psychometric::alpha.CI