## ----setup, include=FALSE-----------------------------------------------------
library(OrdinalCompositions)

data("CitySatisf", package = "OrdinalCompositions")

satisf <- as.matrix(CitySatisf[,14:17]); ya <- satisf/rowSums(satisf)
safety <- as.matrix(CitySatisf[,22:25]); x1a <- safety/rowSums(safety)
air <- as.matrix(CitySatisf[,10:13]); x2a <- air/rowSums(air)

N <- dim(ya)[1]
Cy <- dim(ya)[2]
Cx1 <- dim(x1a)[2]
Cx2 <- dim(x2a)[2]

y <- matrix(0, nrow=N, ncol=Cy)
x1 <- matrix(0, nrow=N, ncol=Cx1)
x2 <- matrix(0, nrow=N, ncol=Cx2)

for(i in 1:N){
  y[i,]  <- ya[i,4:1]
  x1[i,] <- x1a[i,4:1]
  x2[i,] <- x2a[i,4:1]
}

xdata <- list()
ydata <- list()

for(i in 1:N){
  xdata[[i]] <- list(x1[i,], x2[i,])
  ydata[[i]] <- y[i,]
}

a <- b <- weights_orig <- c(2,1,2)

## ----model--------------------------------------------------------------------
res <- ordinal_regression_simplex(
  xdata,
  ydata,
  weights = weights_orig,
  weights_product = list(a,b),
  return_tensor_index = TRUE,
  lambda1_grid = 0,
  compute_opi = TRUE,
  compute_R2 = TRUE,
  compute_OCC = TRUE
)

## ----results------------------------------------------------------------------
cat("Estimated A_hat\n")
A<-res$A_hat
colnames(A) <- apply(res$tensor_index, 1, function(v) paste0("(", v[1], ",", v[2], ")"))
print(round(A, 4))

cat("Wasserstein R2\n")
print(res$R2)

cat("Ordinal Preservation Index (OPI)\n")
print(res$OPI)

cat("Ordinal Compositional Correlation (OCC)\n")
print(res$OCC)