####################################################################################
#                    L2 - RELIABILITY AND MEASUREMENT ERROR                        #
#                                    NMST570                                       #
#                             Last update: 04/10/2018                              #
####################################################################################

# Slides available at http://www.cs.cas.cz/martinkova/documents/NMST570_L2slides.pdf
# HW assignment available at http://www.cs.cas.cz/drabinova/documents/NMST570_HW2.pdf

####################################################################################
# Ex 1.2 ###########################################################################
####################################################################################

# for reproducibility
set.seed(123)

# generating true score T ~ N(5, 1), n = 100
true_score <- rnorm(n = 100, mean = 5, sd = 1)
# generating error term e ~ N(0, 1), n = 100
error <- rnorm(n = 100, mean = 0, sd = 1)

# calculation observed score, X = T + e
observed_score <- true_score + error

# estimating reliability
### formula 1: rel(X) = cor(T, X)^2
(rel1 <- cor(true_score, observed_score)^2)
### formula 2: rel(X) = var(T)/var(X)
(rel2 <- var(true_score)/var(observed_score))
# they are different due to the sampling

# the same, with n = 1000
# generating true score T ~ N(5, 1), n = 100
true_score <- rnorm(n = 1000, mean = 5, sd = 1)
# generating error term e ~ N(0, 1), n = 100
error <- rnorm(n = 1000, mean = 0, sd = 1)

# calculation observed score, X = T + e
observed_score <- true_score + error

# estimating reliability
### formula 1: rel(X) = cor(T, X)^2
(rel1 <- cor(true_score, observed_score)^2)
### formula 2: rel(X) = var(T)/var(X)
(rel2 <- var(true_score)/var(observed_score))
# with increasing sample size, the reliability estimates should be closer to 0.5 (theoretical value)

####################################################################################
# Ex 1.3 ###########################################################################
####################################################################################

library(psychometric)

# reliability of original test
relX <- 0.7
# number of items in original test
I <- 5

#-----------------------------------
# number of items in new test
Istar1 <- 10
m1 <- Istar1/I

# calculating reliability of the new test
(relXstar1 <- SBrel(Nlength = m1, rxx = relX))

#-----------------------------------
# desired reliability
relXstar2 <- 0.9

(m2 <- SBlength(rxxp = relXstar2, rxx = relX))
# calculating number of items in new test
ceiling(Istar2 <- m2*I)