---
title: "Ridge model"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Ridge model 2}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r}
library(rvfl)
```

# Example 1: MPG Prediction (mtcars dataset)

## Load and prepare data
```{r}
data(mtcars)

set.seed(1243)
train_idx <- sample(nrow(mtcars), size = floor(0.8 * nrow(mtcars)))
train_data <- mtcars[train_idx, ]
test_data <- mtcars[-train_idx, -1]
```

## Fit models
```{r}
# Fit regular linear model
start <- proc.time()[3]
lm_model <- lm(mpg ~ ., data = train_data)
print(proc.time()[3] - start)
print(summary(lm_model))
print(confint(lm_model))

# Fit calibrated model 
start <- proc.time()[3]
ridge_model <- rvfl::calibmodel(lambda=10**seq(-10, 10, length.out=100), x = as.matrix(train_data[,-1]), y = train_data$mpg)
print(proc.time()[3] - start)
print(summary(ridge_model))
```

## Make predictions
```{r eval=TRUE}
lm_pred <- predict(lm_model, newdata = test_data, interval = "prediction")
ridge_pred <- predict(ridge_model, newdata = as.matrix(test_data), method = "gaussian")
```

## Compare predictions
```{r eval=TRUE, fig.width=7.5}
results <- data.frame(
  Actual = mtcars[-train_idx, ]$mpg,
  LM_Pred = lm_pred[,"fit"],
  LM_Lower = lm_pred[,"lwr"],
  LM_Upper = lm_pred[,"upr"],
  Ridge_Pred = ridge_pred[,"fit"],
  Ridge_Lower = ridge_pred[,"lwr"], 
  Ridge_Upper = ridge_pred[,"upr"]
)

# Print results
print("Prediction Intervals Comparison:")
print(head(results))

# Calculate coverage and Winkler scores
lm_coverage <- mean(mtcars[-train_idx, ]$mpg >= results$LM_Lower & 
                   mtcars[-train_idx, ]$mpg <= results$LM_Upper)
ridge_coverage <- mean(mtcars[-train_idx, ]$mpg >= results$Ridge_Lower & 
                      mtcars[-train_idx, ]$mpg <= results$Ridge_Upper)

lm_winkler <- misc::winkler_score(mtcars[-train_idx, ]$mpg, results$LM_Lower, results$LM_Upper)
ridge_winkler <- misc::winkler_score(mtcars[-train_idx, ]$mpg, results$Ridge_Lower, results$Ridge_Upper)

print(sprintf("\nPrediction interval metrics:"))
print(sprintf("Linear Model: %.1f%% coverage, %.3f Winkler score", 
              100 * lm_coverage, mean(lm_winkler)))
print(sprintf("Calibrated Model: %.1f%% coverage, %.3f Winkler score", 
              100 * ridge_coverage, mean(ridge_winkler)))

# Set common y-axis limits for both plots
y_limits <- range(c(results$LM_Lower, results$LM_Upper,
                   results$Ridge_Lower, results$Ridge_Upper))

# Plot prediction intervals
par(mfrow=c(1,2))

# Linear Model Plot
plot(results$Actual, results$LM_Pred, 
     main="Linear Model Predictions",
     xlab="Actual MPG", ylab="Predicted MPG",
     ylim=y_limits)
# Add shaded prediction intervals
x_ordered <- order(results$Actual)
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
        c(results$LM_Lower[x_ordered], rev(results$LM_Upper[x_ordered])),
        col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$LM_Pred)  # Replot points over shading
abline(0, 1, col="red", lty=2)  # Add diagonal line

# Ridge Model Plot
plot(results$Actual, results$Ridge_Pred,
     main="Ridge Model Predictions",
     xlab="Actual MPG", ylab="Predicted MPG",
     ylim=y_limits)
# Add shaded prediction intervals
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
        c(results$Ridge_Lower[x_ordered], rev(results$Ridge_Upper[x_ordered])),
        col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$Ridge_Pred)  # Replot points over shading
abline(0, 1, col="red", lty=2)  # Add diagonal line
```

```{r, fig.width=7.5}
# Add simulation plot
par(mfrow=c(1,1))
# Generate 100 simulations
sims <- simulate(ridge_model, newdata = as.matrix(test_data), nsim = 500)
# Plot simulations
matplot(sims, type = "l", 
        col = rgb(0, 0, 1, 0.1), lty = 1,
        xlab = "obs. #", ylab = "Simulated MPG",
        main = "Ridge Model Simulations")
lines(mtcars[-train_idx, ]$mpg, col = "red")        
```

# Example 2: Boston Housing Price Prediction

## Load and prepare data
```{r}
library(MASS)
data(Boston)

set.seed(1243)
train_idx <- sample(nrow(Boston), size = floor(0.8 * nrow(Boston)))
train_data <- Boston[train_idx, ]
test_data <- Boston[-train_idx, -14]  # -14 removes 'medv' (target variable)
```

## Fit models
```{r}
# Fit regular linear model
start <- proc.time()[3]
lm_model <- lm(medv ~ ., data = train_data)
print(proc.time()[3] - start)
print(summary(lm_model))
print(confint(lm_model))

# Fit calibrated model 
start <- proc.time()[3]
ridge_model <- rvfl::calibmodel(lambda=10**seq(-10, 10, length.out=100), x = as.matrix(train_data[,-14]), y = train_data$medv)
print(proc.time()[3] - start)
print(summary(ridge_model$model))
print(confint(ridge_model$model))
#print(simulate(ridge_model, newdata = test_data))
```

## Make predictions and compare
```{r eval=TRUE, fig.width=7.5}
lm_pred <- predict(lm_model, newdata = test_data, interval = "prediction")
ridge_pred <- predict(ridge_model, newdata = as.matrix(test_data), method = "gaussian")

results <- data.frame(
  Actual = Boston[-train_idx, ]$medv,
  LM_Pred = lm_pred[,"fit"],
  LM_Lower = lm_pred[,"lwr"],
  LM_Upper = lm_pred[,"upr"],
  Ridge_Pred = ridge_pred[,"fit"],
  Ridge_Lower = ridge_pred[,"lwr"], 
  Ridge_Upper = ridge_pred[,"upr"]
)

# Print results
print("Prediction Intervals Comparison:")
print(head(results))

# Calculate coverage and Winkler scores
lm_coverage <- mean(Boston[-train_idx, ]$medv >= results$LM_Lower & 
                   Boston[-train_idx, ]$medv <= results$LM_Upper)
ridge_coverage <- mean(Boston[-train_idx, ]$medv >= results$Ridge_Lower & 
                      Boston[-train_idx, ]$medv <= results$Ridge_Upper)

lm_winkler <- misc::winkler_score(Boston[-train_idx, ]$medv, results$LM_Lower, results$LM_Upper)
ridge_winkler <- misc::winkler_score(Boston[-train_idx, ]$medv, results$Ridge_Lower, results$Ridge_Upper)

print(sprintf("\nPrediction interval metrics:"))
print(sprintf("Linear Model: %.1f%% coverage, %.3f Winkler score", 
              100 * lm_coverage, mean(lm_winkler)))
print(sprintf("Calibrated Model: %.1f%% coverage, %.3f Winkler score", 
              100 * ridge_coverage, mean(ridge_winkler)))

# Visualization
# Set common y-axis limits for both plots
y_limits <- range(c(results$LM_Lower, results$LM_Upper,
                   results$Ridge_Lower, results$Ridge_Upper))

par(mfrow=c(1,2))

# Linear Model Plot
plot(results$Actual, results$LM_Pred, 
     main="Linear Model Predictions",
     xlab="Actual Median Value", ylab="Predicted Median Value",
     ylim=y_limits)
x_ordered <- order(results$Actual)
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
        c(results$LM_Lower[x_ordered], rev(results$LM_Upper[x_ordered])),
        col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$LM_Pred)
abline(0, 1, col="red", lty=2)

# Ridge Model Plot
plot(results$Actual, results$Ridge_Pred,
     main="Ridge Model Predictions",
     xlab="Actual Median Value", ylab="Predicted Median Value",
     ylim=y_limits)
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
        c(results$Ridge_Lower[x_ordered], rev(results$Ridge_Upper[x_ordered])),
        col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$Ridge_Pred)
abline(0, 1, col="red", lty=2)

# Add simulation plot
par(mfrow=c(1,1))
sims <- simulate(ridge_model, newdata = as.matrix(test_data), nsim = 500)
matplot(sims, type = "l", 
        col = rgb(0, 0, 1, 0.1), lty = 1,
        xlab = "obs. #", ylab = "Simulated Median Value",
        main = "Ridge Model Simulations")
lines(Boston[-train_idx, ]$medv, col = "red")
```

# Example 3: Economic Indicators (Longley dataset)

## Load and prepare data
```{r}
data(longley)

set.seed(1243)
train_idx <- sample(nrow(longley), size = floor(0.8 * nrow(longley)))
train_data <- longley[train_idx, ]
test_data <- longley[-train_idx, -7]  # -7 removes 'Employed' (target variable)
```

## Fit models
```{r}
# Fit regular linear model
start <- proc.time()[3]
lm_model <- lm(Employed ~ ., data = train_data)
print(proc.time()[3] - start)
print(summary(lm_model))
print(confint(lm_model))

# Fit calibrated model 
start <- proc.time()[3]
ridge_model <- rvfl::calibmodel(lambda=10**seq(-10, 10, length.out=100), x = as.matrix(train_data[,-7]), y = train_data$Employed)
print(proc.time()[3] - start)
print(summary(ridge_model$model))
print(confint(ridge_model$model))
#print(simulate(ridge_model, newdata = test_data))
```

## Make predictions and compare
```{r eval=TRUE, fig.width=7.5}
lm_pred <- predict(lm_model, newdata = test_data, interval = "prediction")
ridge_pred <- predict(ridge_model, newdata = as.matrix(test_data), method = "gaussian")

results <- data.frame(
  Actual = longley[-train_idx, ]$Employed,
  LM_Pred = lm_pred[,"fit"],
  LM_Lower = lm_pred[,"lwr"],
  LM_Upper = lm_pred[,"upr"],
  Ridge_Pred = ridge_pred[,"fit"],
  Ridge_Lower = ridge_pred[,"lwr"], 
  Ridge_Upper = ridge_pred[,"upr"]
)

# Print results
print("Prediction Intervals Comparison:")
print(head(results))

# Calculate coverage and Winkler scores
lm_coverage <- mean(longley[-train_idx, ]$Employed >= results$LM_Lower & 
                   longley[-train_idx, ]$Employed <= results$LM_Upper)
ridge_coverage <- mean(longley[-train_idx, ]$Employed >= results$Ridge_Lower & 
                      longley[-train_idx, ]$Employed <= results$Ridge_Upper)

lm_winkler <- misc::winkler_score(longley[-train_idx, ]$Employed, results$LM_Lower, results$LM_Upper)
ridge_winkler <- misc::winkler_score(longley[-train_idx, ]$Employed, results$Ridge_Lower, results$Ridge_Upper)

print(sprintf("\nPrediction interval metrics:"))
print(sprintf("Linear Model: %.1f%% coverage, %.3f Winkler score", 
              100 * lm_coverage, mean(lm_winkler)))
print(sprintf("Calibrated Model: %.1f%% coverage, %.3f Winkler score", 
              100 * ridge_coverage, mean(ridge_winkler)))

# Visualization
# Set common y-axis limits for both plots
y_limits <- range(c(results$LM_Lower, results$LM_Upper,
                   results$Ridge_Lower, results$Ridge_Upper))

par(mfrow=c(1,2))

# Linear Model Plot
plot(results$Actual, results$LM_Pred, 
     main="Linear Model Predictions",
     xlab="Actual Employment", ylab="Predicted Employment",
     ylim=y_limits)
x_ordered <- order(results$Actual)
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
        c(results$LM_Lower[x_ordered], rev(results$LM_Upper[x_ordered])),
        col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$LM_Pred)
abline(0, 1, col="red", lty=2)

# Ridge Model Plot
plot(results$Actual, results$Ridge_Pred,
     main="Ridge Model Predictions",
     xlab="Actual Employment", ylab="Predicted Employment",
     ylim=y_limits)
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
        c(results$Ridge_Lower[x_ordered], rev(results$Ridge_Upper[x_ordered])),
        col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$Ridge_Pred)
abline(0, 1, col="red", lty=2)

# Add simulation plot
par(mfrow=c(1,1))
sims <- simulate(ridge_model, newdata = as.matrix(test_data), nsim = 500)
matplot(sims, type = "l", 
        col = rgb(0, 0, 1, 0.1), lty = 1,
        xlab = "obs. #", ylab = "Simulated Employment",
        main = "Ridge Model Simulations")
lines(longley[-train_idx, ]$Employed, col = "red")
```

# Example 4: US Crime Rate Analysis

## Load and prepare data
```{r}
data(UScrime)

set.seed(1243)
train_idx <- sample(nrow(UScrime), size = floor(0.8 * nrow(UScrime)))
train_data <- UScrime[train_idx, ]
test_data <- UScrime[-train_idx, -16]  # -16 removes 'y' (crime rate)
```

## Fit models
```{r}
# Fit regular linear model
start <- proc.time()[3]
lm_model <- lm(y ~ ., data = train_data)
print(proc.time()[3] - start)
print(summary(lm_model))
print(confint(lm_model))

# Fit calibrated model 
start <- proc.time()[3]
ridge_model <- rvfl::calibmodel(lambda=10**seq(-10, 10, length.out=100), x = as.matrix(train_data[,-16]), y = train_data$y)
print(proc.time()[3] - start)
print(summary(ridge_model$model))
print(confint(ridge_model$model))
#print(simulate(ridge_model, newdata = test_data))
```

## Make predictions and compare
```{r eval=TRUE, fig.width=7.5}
lm_pred <- predict(lm_model, newdata = test_data, interval = "prediction")
ridge_pred <- predict(ridge_model, newdata = as.matrix(test_data), method = "gaussian")

results <- data.frame(
  Actual = UScrime[-train_idx, ]$y,
  LM_Pred = lm_pred[,"fit"],
  LM_Lower = lm_pred[,"lwr"],
  LM_Upper = lm_pred[,"upr"],
  Ridge_Pred = ridge_pred[,"fit"],
  Ridge_Lower = ridge_pred[,"lwr"], 
  Ridge_Upper = ridge_pred[,"upr"]
)

# Print results
print("Prediction Intervals Comparison:")
print(results)

# Calculate coverage and Winkler scores
lm_coverage <- mean(UScrime[-train_idx, ]$y >= results$LM_Lower & 
                   UScrime[-train_idx, ]$y <= results$LM_Upper)
ridge_coverage <- mean(UScrime[-train_idx, ]$y >= results$Ridge_Lower & 
                      UScrime[-train_idx, ]$y <= results$Ridge_Upper)

lm_winkler <- misc::winkler_score(UScrime[-train_idx, ]$y, results$LM_Lower, results$LM_Upper)
ridge_winkler <- misc::winkler_score(UScrime[-train_idx, ]$y, results$Ridge_Lower, results$Ridge_Upper)

print(sprintf("\nPrediction interval metrics:"))
print(sprintf("Linear Model: %.1f%% coverage, %.3f Winkler score", 
              100 * lm_coverage, mean(lm_winkler)))
print(sprintf("Calibrated Model: %.1f%% coverage, %.3f Winkler score", 
              100 * ridge_coverage, mean(ridge_winkler)))

# Visualization
# Set common y-axis limits for both plots
y_limits <- range(c(results$LM_Lower, results$LM_Upper,
                   results$Ridge_Lower, results$Ridge_Upper))

par(mfrow=c(1,2))

# Linear Model Plot
plot(results$Actual, results$LM_Pred, 
     main="Linear Model Predictions",
     xlab="Actual Crime Rate", ylab="Predicted Crime Rate",
     ylim=y_limits)
x_ordered <- order(results$Actual)
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
        c(results$LM_Lower[x_ordered], rev(results$LM_Upper[x_ordered])),
        col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$LM_Pred)
abline(0, 1, col="red", lty=2)

# Ridge Model Plot
plot(results$Actual, results$Ridge_Pred,
     main="Ridge Model Predictions",
     xlab="Actual Crime Rate", ylab="Predicted Crime Rate",
     ylim=y_limits)
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
        c(results$Ridge_Lower[x_ordered], rev(results$Ridge_Upper[x_ordered])),
        col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$Ridge_Pred)
abline(0, 1, col="red", lty=2)

# Add simulation plot
par(mfrow=c(1,1))
sims <- simulate(ridge_model, newdata = as.matrix(test_data), nsim = 500)
matplot(sims, type = "l", 
        col = rgb(0, 0, 1, 0.1), lty = 1,
        xlab = "obs. #", ylab = "Simulated Crime Rate",
        main = "Ridge Model Simulations")
lines(UScrime[-train_idx, ]$y, col = "red")
```

# Example 5: Car Price Analysis (Cars93 dataset)

## Load and prepare data
```{r}
data(Cars93, package = "MASS")

# Remove rows with missing values
Cars93 <- na.omit(Cars93)

# Select numeric predictors and price as response
predictors <- c("MPG.city", "MPG.highway", "EngineSize", "Horsepower", 
                "RPM", "Rev.per.mile", "Fuel.tank.capacity", "Length", 
                "Wheelbase", "Width", "Turn.circle", "Weight")
car_data <- Cars93[, c(predictors, "Price")]

set.seed(1243)
train_idx <- sample(nrow(car_data), size = floor(0.8 * nrow(car_data)))
train_data <- car_data[train_idx, ]
test_data <- car_data[-train_idx, -which(names(car_data) == "Price")]
```

## Fit models
```{r}
# Fit regular linear model
start <- proc.time()[3]
lm_model <- lm(Price ~ ., data = train_data)
print(proc.time()[3] - start)
print(summary(lm_model))
print(confint(lm_model))

# Fit calibrated model 
start <- proc.time()[3]
ridge_model <- rvfl::calibmodel(lambda=10**seq(-10, 10, length.out=100), x = as.matrix(train_data[,-which(names(train_data) == "Price")]), 
                               y = train_data$Price)
print(proc.time()[3] - start)
print(summary(ridge_model$model))
print(confint(ridge_model$model))
#print(simulate(ridge_model, newdata = as.matrix(test_data)))
```

## Make predictions and compare
```{r eval=TRUE, fig.width=7.5}
lm_pred <- predict(lm_model, newdata = test_data, interval = "prediction")
ridge_pred <- predict(ridge_model, newdata = as.matrix(test_data), method = "gaussian")

results <- data.frame(
  Actual = car_data[-train_idx, "Price"],
  LM_Pred = lm_pred[,"fit"],
  LM_Lower = lm_pred[,"lwr"],
  LM_Upper = lm_pred[,"upr"],
  Ridge_Pred = ridge_pred[,"fit"],
  Ridge_Lower = ridge_pred[,"lwr"], 
  Ridge_Upper = ridge_pred[,"upr"]
)

# Print results
print("Prediction Intervals Comparison:")
print(results)

# Calculate coverage and Winkler scores
lm_coverage <- mean(car_data[-train_idx, "Price"] >= results$LM_Lower & 
                   car_data[-train_idx, "Price"] <= results$LM_Upper)
ridge_coverage <- mean(car_data[-train_idx, "Price"] >= results$Ridge_Lower & 
                      car_data[-train_idx, "Price"] <= results$Ridge_Upper)

lm_winkler <- misc::winkler_score(car_data[-train_idx, "Price"], results$LM_Lower, results$LM_Upper)
ridge_winkler <- misc::winkler_score(car_data[-train_idx, "Price"], results$Ridge_Lower, results$Ridge_Upper)

print(sprintf("\nPrediction interval metrics:"))
print(sprintf("Linear Model: %.1f%% coverage, %.3f Winkler score", 
              100 * lm_coverage, mean(lm_winkler)))
print(sprintf("Calibrated Model: %.1f%% coverage, %.3f Winkler score", 
              100 * ridge_coverage, mean(ridge_winkler)))

# Visualization
# Set common y-axis limits for both plots
y_limits <- range(c(results$LM_Lower, results$LM_Upper,
                   results$Ridge_Lower, results$Ridge_Upper))

par(mfrow=c(1,2))

# Linear Model Plot
plot(results$Actual, results$LM_Pred, 
     main="Linear Model Predictions",
     xlab="Actual Price ($1000s)", ylab="Predicted Price ($1000s)",
     ylim=y_limits)
x_ordered <- order(results$Actual)
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
        c(results$LM_Lower[x_ordered], rev(results$LM_Upper[x_ordered])),
        col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$LM_Pred)
abline(0, 1, col="red", lty=2)

# Ridge Model Plot
plot(results$Actual, results$Ridge_Pred,
     main="Ridge Model Predictions",
     xlab="Actual Price ($1000s)", ylab="Predicted Price ($1000s)",
     ylim=y_limits)
polygon(c(results$Actual[x_ordered], rev(results$Actual[x_ordered])),
        c(results$Ridge_Lower[x_ordered], rev(results$Ridge_Upper[x_ordered])),
        col=rgb(0, 0, 1, 0.2), border=NA)
points(results$Actual, results$Ridge_Pred)
abline(0, 1, col="red", lty=2)

# Add simulation plot
par(mfrow=c(1,1))
sims <- simulate(ridge_model, newdata = as.matrix(test_data), nsim = 500)
matplot(sims, type = "l", 
        col = rgb(0, 0, 1, 0.1), lty = 1,
        xlab = "obs. #", ylab = "Simulated Price ($1000s)",
        main = "Ridge Model Simulations")
lines(car_data[-train_idx, "Price"], col = "red")
```