## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")

## ----fig.width=5, fig.height=5------------------------------------------------
library("Gifi")
library("MPsychoR")
data("granularity")
granularity1 <- scale(granularity[,1:2]) |> as.data.frame()
head(granularity1)
plot(granularity1[,2:1], main = "Scatterplot")

## -----------------------------------------------------------------------------
fitlin1 <- lm(gran ~ -1 + age, data = granularity1)
coef(fitlin1)

## -----------------------------------------------------------------------------
xknots_age <- knotsGifi(granularity$age, type = "E")      
yknots_gran <- knotsGifi(granularity$gran, type = "E")
fitlin2 <-  morals(x = granularity$age, y = granularity$gran, 
   xknots = xknots_age, yknots = yknots_gran, xdegrees = 1, ydegrees = 1,
   xordinal = FALSE, yordinal = FALSE)
fitlin2

## -----------------------------------------------------------------------------
fitquad1 <- lm(gran ~ age + I(age^2), data = granularity)
fitquad1
fitquad2 <-  morals(x = granularity$age, y = granularity$gran, 
  xknots = xknots_age, yknots = yknots_gran, xdegrees = 2, ydegrees = 1, 
  xordinal = FALSE, yordinal = FALSE) 
fitquad2

## ----echo = 2:6, fig.width=5, fig.height=7------------------------------------
op <- par(mfrow = c(2,1))
plot(fitquad2$xhat, fitquad2$yhat, xlab = "Age (transformed)", ylab = "Granularity (transformed)", 
     main = "Optimally Scaled Scatterplot")
lines(fitquad2$xhat, fitquad2$ypred, col = "coral4", lwd = 2)
plot(granularity$age, fitquad2$yhat, xlab = "Age", ylab = "Granularity (transformed)", 
     main = "Quadratic Morals Fit")
ind <- order(granularity$age)
lines(granularity$age[ind], fitquad2$ypred[ind], col = "coral4", lwd = 2)
par(op)

## -----------------------------------------------------------------------------
xknots_age2 <- knotsGifi(granularity$age, "Q", n = 1)
xknots_age2

## -----------------------------------------------------------------------------
xknots_age2[[1]] <- 18

## -----------------------------------------------------------------------------
fitpiece <-  morals(x = granularity$age, y = granularity$gran, 
   xknots = xknots_age2, yknots = yknots_gran, xdegrees = 1, ydegrees = 1,
   xordinal = FALSE, yordinal = FALSE)
fitpiece

## ----echo = 2:6, fig.width=5, fig.height=7------------------------------------
op <- par(mfrow = c(2,1))
plot(fitpiece$xhat, fitpiece$yhat, xlab = "Age (transformed)", ylab = "Granularity (transformed)", 
     main = "Optimally Scaled Scatterplot")
lines(fitpiece$xhat, fitpiece$ypred, col = "coral4", lwd = 2)
plot(granularity$age, fitpiece$yhat, xlab = "Age", ylab = "Granularity (transformed)", 
     main = "Piecewise Linear Morals Fit")
ind <- order(granularity$age)
lines(granularity$age[ind], fitpiece$ypred[ind], col = "coral4", lwd = 2)
par(op)

## -----------------------------------------------------------------------------
xknots_age3 <- knotsGifi(granularity$age, "Q", n = 3)
xknots_age3
fitspline <-  morals(granularity$age, granularity$gran, 
   xknots = xknots_age3, yknots = yknots_gran, xdegrees = 2, ydegrees = 1,
   xordinal = FALSE, yordinal = FALSE)
fitspline

## ----echo = 2:6, fig.width=5, fig.height=7------------------------------------
op <- par(mfrow = c(2,1))
plot(fitspline$xhat, fitspline$yhat, xlab = "Age (transformed)", ylab = "Granularity (transformed)", 
     main = "Optimally Scaled Scatterplot")
lines(fitspline$xhat, fitspline$ypred, col = "coral4", lwd = 2)
plot(granularity$age, fitspline$yhat, xlab = "Age", ylab = "Granularity (transformed)", 
     main = "Spline Morals Fit")
ind <- order(granularity$age)
lines(granularity$age[ind], fitspline$ypred[ind], col = "coral4", lwd = 2)
par(op)

## -----------------------------------------------------------------------------
xknots_age4 <- knotsGifi(granularity$age, "D")
fitmono <-  morals(x = granularity$age, y = granularity$gran, 
  xknots = xknots_age4, yknots = yknots_gran, ydegrees = 1, yordinal = FALSE)

## ----echo = 2:6, fig.width=5, fig.height=7------------------------------------
op <- par(mfrow = c(2,1))
plot(fitmono$xhat, fitmono$yhat, xlab = "Age (transformed)", ylab = "Granularity (transformed)", 
     main = "Optimally Scaled Scatterplot")
lines(fitmono$xhat, fitmono$ypred, col = "coral4", lwd = 2)
plot(granularity$age, fitmono$yhat, xlab = "Age", ylab = "Granularity (transformed)", 
     main = "Monotone Morals Fit")
sfun <- stepfun(sort(granularity$age)[-nrow(granularity)], fitmono$ypred[order(granularity$age)])
plot(sfun, col = "coral4", add = TRUE, pch = 19, cex = 0.7, lwd = 2)
par(op)

## -----------------------------------------------------------------------------
xknots_age5 <- knotsGifi(granularity$age, "D")
fitnom <-  morals(granularity$age, granularity$gran, 
   xknots = xknots_age5, yknots = yknots_gran, xdegrees = 1, ydegrees = 1,
   xordinal = FALSE, yordinal = FALSE)

## ----echo = 2:6, fig.width=5, fig.height=7------------------------------------
op <- par(mfrow = c(2,1))
plot(fitnom$xhat, fitnom$yhat, xlab = "Age (transformed)", ylab = "Granularity (transformed)", 
     main = "Optimally Scaled Scatterplot")
lines(fitnom$xhat, fitnom$ypred, col = "coral4", lwd = 2)
plot(granularity$age, fitnom$yhat, xlab = "Age", ylab = "Granularity (transformed)", 
     main = "Nominal Morals Fit")
ind <- order(granularity$age)
lines(granularity$age[ind], fitnom$ypred[ind], col = "coral4", lwd = 2)
par(op)

## ----cvmorals, cache=TRUE-----------------------------------------------------
set.seed(123)
cvlin <- cv(fitlin2, folds = 10)
cvquad <- cv(fitquad2, folds = 10)
cvpiece <- cv(fitpiece, folds = 10)
cvspline <- cv(fitspline, folds = 10)
cvmono <- cv(fitmono, folds = 10)
cvnom <- cv(fitnom, folds = 10)
cvvec <- c(cvlin, cvquad, cvpiece, cvspline, cvmono, cvnom)
r2vec <- c(fitlin2$smc, fitquad2$smc, fitpiece$smc, fitspline$smc,
          fitmono$smc, fitnom$smc)
cvr2 <- cbind(cvvec, r2vec)
dimnames(cvr2) <- list(c("linear", "quadratic", "piecewise", "spline",
                         "monotone", "nominal"), c("CV-error", "R2"))
round(cvr2, 5)

## -----------------------------------------------------------------------------
granularity2 <- granularity 
granularity2$gender <- as.numeric(granularity$gender)-1
granularity2 <- scale(granularity2) |> as.data.frame()
head(granularity2)
fitmlin1 <- lm(gran ~ -1 + age*gender, data = granularity2)
fitmlin1

## -----------------------------------------------------------------------------
granularity2$int <- granularity2$age * granularity2$gender
xknots_age <- knotsGifi(granularity2[,2:4], "E")  
yknots_gran <- knotsGifi(granularity2$gran, "E")
fitmlin2 <-  morals(x = granularity2[,2:4], y= granularity2$gran, 
   xknots = xknots_age, yknots = yknots_gran, xdegrees = 1, ydegrees = 1,
   xordinal = FALSE, yordinal = FALSE)
fitmlin2

## -----------------------------------------------------------------------------
library("MASS")
grancat <- cut(granularity$gran, 5, labels = 1:5)
fitord1 <- polr(grancat ~ age + I(age^2) + gender, data = granularity)
summary(fitord1)

## -----------------------------------------------------------------------------
granularity3 <- granularity 
granularity3$gender <- as.numeric(granularity$gender)
xknots_age <- knotsGifi(granularity3[,2:3], "E")  
yknots_gran2 <- knotsGifi(grancat, "D")
fitord2 <-  morals(x = granularity3[,2:3], y = as.numeric(grancat),
  xknots = xknots_age, yknots = yknots_gran2, xdegrees = c(2, -1), ydegrees = 1, 
  xordinal = FALSE, yordinal = TRUE) 
fitord2

## ----fig.height=7, fig.width=7------------------------------------------------
plot(fitord2, "transplot", main = c("Granularity Categorical", "Age", "Gender"))

## ----echo = 2:11, fig.width=5, fig.height=7-----------------------------------
op <- par(mfrow = c(2,1))
plot(fitord2$xhat[,1], fitord2$yhat, xlab = "Age (transformed)", ylab = "Granularity (transformed)", 
     main = "Optimally Scaled Scatterplot")
lines(fitord2$xhat[,1][granularity3$gender == 1], fitord2$ypred[granularity3$gender == 1], col = "coral4", lwd = 2)
lines(fitord2$xhat[,1][granularity3$gender == 2], fitord2$ypred[granularity3$gender == 2], col = "cadetblue", lwd = 2)
legend(-0.02, 0.14, bty = "n", legend = c("male", "female"), lty = 1, col = c("coral4", "cadetblue"))
plot(granularity3$age, fitord2$yhat, xlab = "Age", ylab = "Granularity (transformed)", main = "Ordinal Polynomial Morals Fit")
ind1 <- order(granularity3$age[granularity3$gender == 1])
lines(granularity3$age[granularity3$gender == 1][ind1], fitord2$ypred[granularity3$gender == 1][ind1], col = "coral4", lwd = 2)
ind2 <- order(granularity3$age[granularity3$gender == 2])
lines(granularity3$age[granularity3$gender == 2][ind2], fitord2$ypred[granularity3$gender == 2][ind2], col = "cadetblue", lwd = 2)
legend(20, 0.14, bty = "n", legend = c("male", "female"), lty = 1, col = c("coral4", "cadetblue"))
par(op)