## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----message=FALSE------------------------------------------------------------
library(mvord)

## -----------------------------------------------------------------------------
data("data_cr", package = "mvord")
str(data_cr, vec.len = 3)
head(data_cr, n = 3)

## ----eval = FALSE-------------------------------------------------------------
# res_cor_probit_2raters <- mvord(formula = MMO2(rater1, rater2) ~ 0 + LR + LEV + PR + RSIZE + BETA,
#                                threshold.constraints = c(1, 1),
#                                data = data_cr)

## ----echo = FALSE, results = 'hide', eval = TRUE------------------------------
cache <- TRUE
FILE <- "res_cor_probit_2raters.rda"
if (cache & file.exists(FILE)) {
  load(FILE)
} else {
  if (cache) {
res_cor_probit_2raters <- mvord(formula = MMO2(rater1, rater2) ~ 0 + LR + LEV + PR + RSIZE + BETA,
                               threshold.constraints = c(1, 1),
                               data = data_cr)
res_cor_probit_2raters <- mvord:::reduce_size.mvord(res_cor_probit_2raters)
  save(res_cor_probit_2raters, file  = FILE)
  } else {
      if(file.exists(FILE)) file.remove(FILE)
  }
}

## -----------------------------------------------------------------------------
summary(res_cor_probit_2raters)

## -----------------------------------------------------------------------------
coef(res_cor_probit_2raters)

## ----eval = FALSE-------------------------------------------------------------
# res_cor_logit_3raters <- mvord(formula = MMO2(rater1, rater2, rater3) ~ 0 +
#                                  LR + LEV + PR + RSIZE + BETA,
#                                coef.constraints = c(1, 1, 1),
#                                data = data_cr,
#                                link = mvlogit())

## ----echo = FALSE, results = 'hide', eval = TRUE------------------------------
cache <- TRUE
FILE <- "res_cor_logit_3raters.rda"
if (cache & file.exists(FILE)) {
  load(FILE)
} else {
  if (cache) {
res_cor_logit_3raters <- mvord(formula = MMO2(rater1, rater2, rater3) ~ 0 + LR + LEV + PR + RSIZE + BETA,
                               coef.constraints = c(1, 1, 1),
                               data = data_cr,
                               link = mvlogit())
res_cor_logit_3raters <- mvord:::reduce_size.mvord(res_cor_logit_3raters)
  save(res_cor_logit_3raters, file  = FILE)
  } else {
      if(file.exists(FILE)) file.remove(FILE)
  }
}

## -----------------------------------------------------------------------------
summary(res_cor_logit_3raters)

## -----------------------------------------------------------------------------
thresholds(res_cor_logit_3raters)

## ----eval = FALSE-------------------------------------------------------------
# res_cor_probit_simple <- mvord(formula = MMO2(rater1, rater2, rater3, rater4) ~ 0 + LR + LEV + PR + RSIZE + BETA, data = data_cr)

## ----echo = FALSE, results = 'hide', eval = TRUE------------------------------
cache <- TRUE
FILE <- "res_cor_probit_simple.rda"
if (cache & file.exists(FILE)) {
  load(FILE)
} else {
  if (cache) {
res_cor_probit_simple <- mvord(
    formula =
    MMO2(rater1, rater2, rater3, rater4) ~ 0 + LR + LEV + PR + RSIZE + BETA, data = data_cr)
res_cor_probit_simple <- mvord:::reduce_size.mvord(res_cor_probit_simple)
  save(res_cor_probit_simple, file  = FILE)
  } else {
      if(file.exists(FILE)) file.remove(FILE)
  }
}

## -----------------------------------------------------------------------------
summary(res_cor_probit_simple)

## -----------------------------------------------------------------------------
thresholds(res_cor_probit_simple)

## -----------------------------------------------------------------------------
coef(res_cor_probit_simple)

## -----------------------------------------------------------------------------
error_structure(res_cor_probit_simple)[[11]]

## ----eval = FALSE-------------------------------------------------------------
# res_cor_logit <- mvord(formula = MMO2(rater1, rater2, rater3, rater4) ~
#     0 + LR + LEV + PR + RSIZE + BETA, data = data_cr, link = mvlogit(),
#     coef.constraints = cbind(LR = c(1, 1, 1, 1),
#                              LEV = c(1, 2, 3, 4),
#                              PR = c(1, 1, 1, 1),
#                              RSIZE = c(1, 1, 1, 2),
#                              BETA = c(1, 1, 2, 3)),
#     threshold.constraints = c(1, 1, 2, 3))

## ----echo = FALSE, results = 'hide', eval = TRUE------------------------------
FILE <- "res_cor_logit.rda"
if (cache & file.exists(FILE)) {
  load(FILE)
} else {
  if (cache) {
res_cor_logit <- mvord(formula = MMO2(rater1, rater2, rater3, rater4) ~
    0 + LR + LEV + PR + RSIZE + BETA, data = data_cr, link = mvlogit(),
  coef.constraints = cbind(
    c(1,1,1,1),
    c(1,2,3,4),
    c(1,1,1,1),
    c(1,1,1,2),
    c(1,1,2,3)),
    threshold.constraints = c(1, 1, 2, 3))
#res_cor_logit <- mvord:::reduce_size2.mvord(res_cor_logit)
  save(res_cor_logit, file  = FILE)
  } else {
      if(file.exists(FILE)) file.remove(FILE)
  }

}

## -----------------------------------------------------------------------------
summary(res_cor_logit)

## -----------------------------------------------------------------------------
constraints(res_cor_logit)

## -----------------------------------------------------------------------------
BIC(res_cor_probit_simple, res_cor_logit)

## -----------------------------------------------------------------------------
AIC(res_cor_probit_simple)
AIC(res_cor_logit)

## -----------------------------------------------------------------------------
logLik(res_cor_probit_simple)
logLik(res_cor_logit)

## -----------------------------------------------------------------------------
mp <- marginal_predict(res_cor_logit, type = "class")

## -----------------------------------------------------------------------------
table(data_cr$rater1, mp$rater1)
table(data_cr$rater2, mp$rater2)
table(data_cr$rater3, mp$rater3)
table(data_cr$rater4, mp$rater4)

## ----eval = TRUE--------------------------------------------------------------
jp <- predict(res_cor_logit, type = "prob")
head(jp)

## ----eval = TRUE--------------------------------------------------------------
jp_AAFL <- joint_probabilities(res_cor_logit, 
                          response.cat = c("A", "A", "F", "L"),  
                          type = "prob")
head(jp_AAFL)

## ----eval = TRUE--------------------------------------------------------------
nc <- 5 # number of classes for rater 1 and 2
jp_mat <- matrix(nrow = nrow(data_cr), ncol = nc)
for (cl in 1:nc) {
  jp_mat[, cl] <- joint_probabilities(
    res_cor_logit, 
    response.cat = c(LETTERS[cl], LETTERS[cl], "F", "L"),  
    type = "prob")
}
jp2 <- rowSums(jp_mat)
head(jp2)

## ----eval = FALSE-------------------------------------------------------------
# jpred <- predict(res_cor_logit, type = "class")

## ----echo = FALSE, results = 'hide', eval = TRUE------------------------------
FILE <- "jp.rda"
if (cache & file.exists(FILE)) {
  load(FILE)
} else {
  if (cache) {
  jpred <- predict(res_cor_logit, type = "class")
  save(jpred, file  = FILE)
  } else {
      if(file.exists(FILE)) file.remove(FILE)
  }
}

## -----------------------------------------------------------------------------
table(data_cr$rater1, jpred$rater1)
table(data_cr$rater2, jpred$rater2)
table(data_cr$rater3, jpred$rater3)
table(data_cr$rater4, jpred$rater4)

## -----------------------------------------------------------------------------
data("data_cr_panel", package = "mvord")
str(data_cr_panel, vec.len = 3)
head(data_cr_panel, n = 3)

## ----eval=F-------------------------------------------------------------------
# res_AR1_probit <- mvord(formula = MMO(rating, firm_id, year) ~ LR + LEV +
#   PR + RSIZE + BETA,
#   error.structure = cor_ar1(~ BSEC), link = mvprobit(),
#   data = data_cr_panel,
#   coef.constraints = c(rep(1, 4), rep(2, 4)),
#   threshold.constraints = rep(1, 8),
#   threshold.values = rep(list(c(0, NA, NA, NA)),8),
#   control = mvord.control(solver = "BFGS"))

## ----echo = FALSE, results = 'hide', eval = TRUE------------------------------
FILE <- "res_AR1_probit.rda"
if (cache & file.exists(FILE)) {
  load(FILE)
} else {
  if (cache) {
res_AR1_probit <- mvord(
  formula = MMO(rating, firm_id, year) ~ LR + LEV + PR + RSIZE +  BETA,
  data = data_cr_panel,
  error.structure = cor_ar1(~ BSEC),
  coef.constraints = c(rep(1,4), rep(2,4)),
  threshold.constraints = c(rep(1,8)),
  threshold.values = rep(list(c(0,NA,NA,NA)),8),
  link = mvprobit(),
  control = mvord.control(solver = "BFGS",  
                          solver.optimx.control = list(trace = TRUE)))
   res_AR1_probit <- mvord:::reduce_size.mvord(res_AR1_probit)
  save(res_AR1_probit, file  = FILE)
  } else {
      if(file.exists(FILE)) file.remove(FILE)
  }
}

## -----------------------------------------------------------------------------
summary(res_AR1_probit, short = TRUE, call = FALSE)

## -----------------------------------------------------------------------------
error_structure(res_AR1_probit)

## -----------------------------------------------------------------------------
head(error_structure(res_AR1_probit, type = "corr"), n = 3)

## -----------------------------------------------------------------------------
head(error_structure(res_AR1_probit, type = "sigmas"), n = 1)

## ----include=FALSE------------------------------------------------------------
# load("data_SRHS_long.rda")

## ----eval=TRUE, include=TRUE--------------------------------------------------
data_path <- system.file("extdata", "data_SRHS_long.rda", package = "mvord")
data_SRHS_long <- get(load(data_path))

## -----------------------------------------------------------------------------
str(data_SRHS_long)

## ----eval=FALSE, include=TRUE-------------------------------------------------
# res_srhs <- mvord(formula = MMO(srhs, id, t) ~ 0 + factor(gender) +
# 	factor(race) + factor(education) + age,
# 	data = data_SRHS_long,
# 	threshold.constraints = rep(1, 8),
# 	coef.constraints = rep(1, 8),
# 	error.structure = cor_ar1(~ 1), link = mvlogit(),
# 	PL.lag = 2)

## ----include=FALSE------------------------------------------------------------
FILE <- "res_srhs.rda"
if (cache & file.exists(FILE)) {
  load(FILE)
} else {
  if (cache) {
    res_srhs <- mvord(formula = MMO(srhs, id, t) ~ 0 +  factor(gender) +
   	    factor(race) + factor(education) + age,
                  data = data_SRHS_long,
                  link = mvlogit(),
                  threshold.constraints = rep(1, 8),
                  coef.constraints = rep(1, 8),
                  error.structure = cor_ar1(~1), PL.lag = 2)
   res_srhs <- mvord:::reduce_size.mvord(res_srhs)
   save(res_srhs, file  = FILE)
  } else {
      if(file.exists(FILE)) file.remove(FILE)
  }
}

## -----------------------------------------------------------------------------
unique(error_structure(res_srhs, type = "corr"))

## -----------------------------------------------------------------------------
summary(res_srhs, call = FALSE)

## ----echo=FALSE, eval=FALSE---------------------------------------------------
# ## links are broken
# N <- "http://www-math.bgsu.edu/~albert/ord_book/Chapter5/essay_data/N.dat"
# X <- "http://www-math.bgsu.edu/~albert/ord_book/Chapter5/essay_data/X.dat"
# y  <- read.delim(url(N), header = F, sep = "")
# wl <- read.delim(url(X), header = F, sep = "")[,2]
# essay_data <- cbind.data.frame(y, wl)
# colnames(essay_data)[1:5] <- paste0("Judge", 1:5)
# save(essay_data, file =  "essay_data.rda")

## ----include=TRUE-------------------------------------------------------------
data("essay_data", package = "mvord")

## -----------------------------------------------------------------------------
head(essay_data)

## ----eval=FALSE, include=TRUE-------------------------------------------------
# res_essay_0 <- mvord(
#   formula = MMO2(Judge1, Judge2, Judge3, Judge4, Judge5) ~ -1,
#   data = essay_data, threshold.constraints = rep(1, 5),
#   coef.constraints = rep(1, 5))

## ----message=FALSE, warning=FALSE, include=FALSE------------------------------
FILE <- "res_essay.rda"
if (cache & file.exists(FILE)) {
  load(FILE)
} else {
  if (cache) {
    res_essay_wl <- mvord(
      formula = MMO2(Judge1, Judge2, Judge3, Judge4, Judge5) ~ 0 + wl,
      data = essay_data, threshold.constraints = rep(1, 5),
      coef.constraints = rep(1, 5))
    res_essay_wl <- mvord:::reduce_size2.mvord(res_essay_wl)
    res_essay_0 <- mvord(
      formula = MMO2(Judge1, Judge2, Judge3, Judge4, Judge5) ~ -1,
      data = essay_data, threshold.constraints = rep(1, 5),
      coef.constraints = rep(1,5))
    res_essay_0 <- mvord:::reduce_size.mvord(res_essay_0)
    save(res_essay_0, res_essay_wl, file  = FILE)
  } else {
      if(file.exists(FILE)) file.remove(FILE)
  }
}

## -----------------------------------------------------------------------------
summary(res_essay_0, call = FALSE)

## ----eval=FALSE, include=TRUE-------------------------------------------------
# res_essay_wl <- mvord(
#   formula = MMO2(Judge1, Judge2, Judge3, Judge4, Judge5) ~ 0 + wl,
#   data = essay_data, threshold.constraints = rep(1, 5),
#   coef.constraints = rep(1, 5))

## -----------------------------------------------------------------------------
summary(res_essay_wl, call = FALSE)

## -----------------------------------------------------------------------------
agree_prob_list <- lapply(1:10, function(i)
  joint_probabilities(res_essay_wl, rep(i, 5)))
agree_prob <- Reduce("+", agree_prob_list)
summary(agree_prob)

## ----fig.align='center', fig.width=6, fig.height=6----------------------------
plot(essay_data$wl, agree_prob,
  xlab = "word length", ylab = "probability of agreement")