####################################################################################
#                 L5 - REGRESSION MODELS FOR ITEM DESCRIPTION                      #
#                                    NMST570                                       #
#                             Last update: 26/10/2018                              #
####################################################################################

# Slides available at http://www.cs.cas.cz/martinkova/NMST570.html
# HW assignment available at http://www.cs.cas.cz/drabinova/documents/NMST570_HW5.pdf

####################################################################################
# Ex. 1.1. #########################################################################
####################################################################################

# Function for logistic model (1)
# x is total score, b0 and b1 are items parameters
# returns probability of correct answer
logistic1 <- function(x, b0, b1){
  exp(b0 + b1*x)/(1 + exp(b0 + b1*x))
}

# Example:
# probability of correct answer for x = 10, b0 = -4, b1 = 0.3 is given by
logistic1(x = 10, b0 = -4, b1 = 0.3)

# Plot of the item for various levels of total score
x <- seq(0, 20, 0.1)                         # levels of total score
prob1 <- logistic1(x = x, b0 = -4, b1 = 0.3) # calculation of probability

plot(prob1 ~ x, type = "l", ylim = c(0, 1))

####################################################################################
# Ex. 1.2. #########################################################################
####################################################################################

# Function for logistic model (2)
# z is standardized total score, b0tilde and b1tilde are parameters
# returns probability of correct answer
logistic2 <- function(z, b0tilde, b1tilde){
  exp(b0tilde + b1tilde*z)/(1 + exp(b0tilde + b1tilde*z))
}

# Example:
# probability of correct answer for z = 1, b0 = 0.5, b1 = 1 is given by
logistic2(z = 1, b0tilde = 0.5, b1tilde = 1)

# Plot of the item for various levels of standardized total score
z <- seq(-4, 4, 0.1)                                  # levels of standardized total score
prob2 <- logistic2(z = z, b0tilde = 0.5, b1tilde = 1) # calculation of probability

plot(prob2 ~ z, type = "l", ylim = c(0, 1))

####################################################################################
# Ex. 1.3. #########################################################################
####################################################################################

# Function for logistic model (3)
# z is standardized total score, a and b are parameters
# returns probability of correct answer
logistic3 <- function(z, a, b){
  exp(a * (z - b))/(1 + exp(a * (z - b)))
}

# Example:
# probability of correct answer for z = -1, a = 1, b = -2 is given by
logistic3(z = -1, a = 1, b = -2)

# Plot of the item for various levels of standardized total score
z <- seq(-4, 4, 0.1)                     # levels of standardized total score
prob3 <- logistic3(z = z, a = 1, b = -2) # calculation of probability

plot(prob3 ~ z, type = "l", ylim = c(0, 1))

####################################################################################
# Ex. 2.1. #########################################################################
####################################################################################

# Function for generalized logistic model (4)
# z is standardized total score, a, b, c and d are parameters
# returns probability of correct answer
logistic4 <- function(z, a, b, c, d){
  c + (d - c) * exp(a * (z - b))/(1 + exp(a * (z - b)))
}

# Example:
# probability of correct answer for z = -1, a = 1, b = -2, c = 0.2, d = 0.95 is given by
logistic4(z = -1, a = 1, b = -2, c = 0.2, d = 0.95)

# Plot of the item for various levels of standardized total score
z <- seq(-4, 4, 0.1)                                        # levels of standardized total score
prob4 <- logistic4(z = z, a = 1, b = -2, c = 0.2, d = 0.95) # calculation of probability

plot(prob4 ~ z, type = "l", ylim = c(0, 1))
abline(a = 0.95, b = 0, col = "red", lty = 2)
abline(a = 0.2, b = 0, col = "red", lty = 2)

####################################################################################
# Ex. 2.2. #########################################################################
####################################################################################

# Function for multinomial model (5)
# z is standardized total score, a and b are vectors of parameters
# returns probability of correct answer
logistic5 <- function(z, a, b){
  k <- length(a)
  res <- matrix(NA, nrow = length(z), ncol = k)
  for (i in 1:k){
    res[, i] <- exp(a[i] * (z - b[i]))
  }
  res <- res/apply(res, 1, sum)
  colnames(res) <- paste0("k = ", 1:k)
  rownames(res) <- paste0("z = ", z)
  t(res)
}

# Example:
# probabilities of option selection for z = 1, a = c(0, 1, 1.5), b = c(0, -1, 0.5) is given by
logistic5(z = 1, a  = c(0, 1, 1.5), b = c(0, -1, 0.5))

# Plot of the item for various levels of standardized total score
z <- seq(-4, 4, 0.1)                                            # levels of standardized total score
prob5 <- logistic5(z = z, a  = c(0, 1, 1.5), b = c(0, -1, 0.5)) # calculation of probabilities

plot(prob5[1, ] ~ z, type = "l", ylim = c(0, 1))
lines(prob5[2, ] ~ z)
abline(v = -1, col = "blue", lty = 2)
lines(prob5[3, ] ~ z)
abline(v = 0.5, col = "blue", lty = 2)
abline(v = 3.5, col = "blue", lty = 2)