title: “Shape Constraints with BsplineQuantReg” author: “Alexandre Abbes” date: “2026-08-20” output: rmarkdown::html_vignette vignette: > % % % —

## Introduction

This vignette covers shape constraints in B-spline quantile regression using the Karlin-Studden SOCP formulation. Shape constraints allow you to incorporate prior knowledge about the function’s behavior (monotonicity, convexity, etc.) into the estimation process.

## Available Demos

The package includes several demos that illustrate shape constraints:

    # Example of available demos
    dev.new()
    demo(package = "BsplineQuantReg",temperature2)

    ##
    ##
    ##  demo(temperature2)
    ##  ---- ~~~~~~~~~~~~
    ##
    ## > # Demo: Temperature Anomaly Trend Analysis
    ## > # This demo shows temperature data with various monotonicity constraints
    ## > # Data source: Global temperature anomalies (1880-1992)
    ## > #
    ## > # Constraint scenarios:
    ## > # 1. Unconstrained fit
    ## > # 2. Full increasing monotonicity (too restrictive)
    ## > # 3. Partial increasing only after 1970
    ## > # 4. Mixed constraints: decreasing (1945-1970) + increasing elsewhere
    ## >
    ## > library(BsplineQuantReg)
    ##
    ## > oldpar <- par(mfrow = c(2,2))
    ##
    ## > cat("========================================\n")
    ## ========================================
    ##
    ## > cat("Demo: Temperature Anomaly Trend Analysis\n")
    ## Demo: Temperature Anomaly Trend Analysis
    ##
    ## > cat("========================================\n\n")
    ## ========================================
    ##
    ##
    ## > if (!exists("degree")){degree=3}
    ##
    ## > # Create temperature data (1880-1992)
    ## > years <- 1880:1992
    ##
    ## > temperature <- c(
    ## +   -0.32, -0.32, -0.40, -0.39, -0.65, -0.43, -0.40, -0.52, -0.30, -0.12,
    ## +   -0.40, -0.42, -0.39, -0.45, -0.35, -0.36, -0.19, -0.14, -0.37, -0.22,
    ## +   0.00, -0.08, -0.24, -0.36, -0.49, -0.27, -0.19, -0.43, -0.29, -0.30,
    ## +   -0.29, -0.29, -0.28, -0.23, -0.04, -0.02, -0.24, -0.42, -0.35, -0.16,
    ## +   -0.17, -0.09, -0.13, -0.16, -0.14, -0.14,  0.10, -0.03,  0.03, -0.18,
    ## +   -0.06,  0.04,  0.02, -0.13,  0.03, -0.06,  0.02,  0.13,  0.13, -0.03,
    ## +   0.15,  0.12,  0.10,  0.04,  0.11, -0.04,  0.01,  0.13, -0.01, -0.06,
    ## +   -0.14, -0.02,  0.04,  0.14, -0.07, -0.06, -0.17,  0.10,  0.10,  0.05,
    ## +   -0.01,  0.08,  0.02,  0.02, -0.26, -0.16, -0.09, -0.02, -0.12,  0.03,
    ## +   0.04, -0.11, -0.07,  0.19, -0.07, -0.05, -0.22,  0.16,  0.09,  0.14,
    ## +   0.28,  0.39,  0.07,  0.29,  0.11,  0.11,  0.16,  0.32,  0.35,  0.25,
    ## +   0.47,  0.41,  0.13
    ## + )
    ##
    ## > cat(sprintf("Years: %d to %d (%d observations)\n",
    ## +             min(years), max(years), length(years)))
    ## Years: 1880 to 1992 (113 observations)
    ##
    ## > cat(sprintf("Temperature range: [%.2f, %.2f] C\n\n",
    ## +             min(temperature), max(temperature)))
    ## Temperature range: [-0.65, 0.47] C
    ##
    ##
    ## > # Normalize years to [0,1] for numerical stability
    ## > xtab <- (years - min(years)) / (max(years) - min(years))
    ##
    ## > ytab <- temperature
    ##
    ## > # Define knots at specific years
    ## > target_knots <- c(1917, 1936, 1945, 1970)
    ##
    ## > knots_years <- sort(c(min(years), target_knots, max(years)))
    ##
    ## > knots <- (knots_years - min(years)) / (max(years) - min(years))
    ##
    ## > cat("Knots at years:", knots_years, "\n\n")
    ## Knots at years: 1880 1917 1936 1945 1970 1992
    ##
    ##
    ## > # Create evaluation grid
    ## > x_eval <- seq(0, 1, length.out = 300)
    ##
    ## > years_eval <- min(years) + x_eval * (max(years) - min(years))
    ##
    ## > # Helper function to get interval index for a given year
    ## > get_interval_idx <- function(year, knots_years) {
    ## +   which(knots_years <= year)[length(which(knots_years <= year))]
    ## + }
    ##
    ## > # ============================================================
    ## > # Model 1: Unconstrained median regression
    ## > # ============================================================
    ## > cat("=== 1. Unconstrained median regression ===\n")
    ## === 1. Unconstrained median regression ===
    ##
    ## > fit_uncon <- quantile_spline(degree=degree,xtab, ytab, knots, tau = 0.5,
    ## +                                     monot = 0, convcons = 0)

    ## Warning in highs::hi_new_solver(model): WARNING: LP matrix packed vector
    ## contains 4 |value| in [1.78251e-11, 1.78572e-11] less than or equal to 1e-09:
    ## ignored

    ##
    ## > y_uncon <- spline_eval(fit_uncon, x_eval)
    ##
    ## > # ============================================================
    ## > # Model 2: Full monotonicity (increasing everywhere)
    ## > # ============================================================
    ## > cat("\n=== 2. Full increasing monotonicity (everywhere) ===\n")
    ##
    ## === 2. Full increasing monotonicity (everywhere) ===
    ##
    ## > fit_full_inc <- quantile_spline(degree=degree,xtab, ytab, knots, tau = 0.5,
    ## +                                        monot = 1, convcons = 0)

    ## Not enough monotonicity constraints, completing with 0

    ## Warning in highs::hi_new_solver(model): WARNING: LP matrix packed vector
    ## contains 9 |value| in [8.88178e-16, 1.78572e-11] less than or equal to 1e-09:
    ## ignored

    ##
    ## > y_full_inc <- spline_eval(fit_full_inc, x_eval)
    ##
    ## > # ============================================================
    ## > # Model 3: Partial monotonicity (only after 1970)
    ## > # ============================================================
    ## > br_year_1970 <- 1970
    ##
    ## > monot_partial <- rep(0, length(knots) - 1)
    ##
    ## > for (i in 1:length(monot_partial)) {
    ## +   interval_end <- knots_years[i + 1]
    ## +   if (interval_end > br_year_1970) {
    ## +     monot_partial[i] <- 1
    ## +   }
    ## + }
    ##
    ## > cat("\n=== 3. Partial monotonicity (increasing only after 1970) ===\n")
    ##
    ## === 3. Partial monotonicity (increasing only after 1970) ===
    ##
    ## > cat("   Constrained intervals:", which(monot_partial == 1), "\n")
    ##    Constrained intervals: 5
    ##
    ## > fit_partial <- quantile_spline(degree=degree,xtab, ytab, knots, tau = 0.5,
    ## +                                       monot = monot_partial, convcons = 0)

    ## Not enough monotonicity constraints, completing with 0

    ## Warning in highs::hi_new_solver(model): WARNING: LP matrix packed vector
    ## contains 5 |value| in [8.88178e-16, 1.78572e-11] less than or equal to 1e-09:
    ## ignored

    ##
    ## > y_partial <- spline_eval(fit_partial, x_eval)
    ##
    ## > # ============================================================
    ## > # Model 4: Mixed constraints
    ## > # Decreasing between 1945 and 1970, increasing elsewhere
    ## > # ============================================================
    ## > br_year_start <- 1945
    ##
    ## > br_year_end <- 1970
    ##
    ## > monot_mixed <- rep(0, length(knots) - 1)
    ##
    ## > for (i in 1:length(monot_mixed)) {
    ## +   interval_start <- knots_years[i]
    ## +   interval_end <- knots_years[i + 1]
    ## +
    ## +   if (interval_start >= br_year_end) {
    ## +     monot_mixed[i] <- 1      # increasing after 1970
    ## +   } else if (interval_end <= br_year_start) {
    ## +     monot_mixed[i] <- 1      # increasing before 1945
    ## +   } else if (interval_start >= br_year_start && interval_end <= br_year_end) {
    ## +     monot_mixed[i] <- -1     # decreasing between 1945 and 1970
    ## +   } else if (interval_start < br_year_start && interval_end > br_year_start) {
    ## +     monot_mixed[i] <- 1      # increasing (partial interval before 1945)
    ## +   } else if (interval_start < br_year_end && interval_end > br_year_end) {
    ## +     monot_mixed[i] <- 1      # increasing (partial interval after 1970)
    ## +   }
    ## + }
    ##
    ## > cat("\n=== 4. Mixed constraints ===\n")
    ##
    ## === 4. Mixed constraints ===
    ##
    ## > cat("   Decreasing between 1945 and 1970\n")
    ##    Decreasing between 1945 and 1970
    ##
    ## > cat("   Increasing elsewhere\n")
    ##    Increasing elsewhere
    ##
    ## > cat("   Constraint vector:", monot_mixed, "\n")
    ##    Constraint vector: 1 1 1 -1 1
    ##
    ## > fit_mixed <- quantile_spline(degree=degree,xtab, ytab, knots, tau = 0.5,
    ## +                                     monot = monot_mixed, convcons = 0)

    ## Not enough monotonicity constraints, completing with 0

    ## Warning in highs::hi_new_solver(model): WARNING: LP matrix packed vector
    ## contains 9 |value| in [8.88178e-16, 1.78572e-11] less than or equal to 1e-09:
    ## ignored

    ##
    ## > y_mixed <- spline_eval(fit_mixed, x_eval)
    ##
    ## > # ============================================================
    ## > # Model 5: Multiple quantiles with mixed constraints
    ## > # ============================================================
    ## > cat("\n=== 5. Multiple quantiles (0.1, 0.5, 0.9) with mixed constraints ===\n")
    ##
    ## === 5. Multiple quantiles (0.1, 0.5, 0.9) with mixed constraints ===
    ##
    ## > tau_multi <- c(0.1, 0.5, 0.9)
    ##
    ## > fits_multi <- list()
    ##
    ## > for (i in seq_along(tau_multi)) {
    ## +   fits_multi[[i]] <- quantile_spline(degree=degree,xtab, ytab, knots,
    ## +                                             tau = tau_multi[i],
    ## +                                             monot = monot_mixed,
    ## +                                             convcons = 0)
    ## +   cat(sprintf("   tau = %.1f done\n", tau_multi[i]))
    ## + }

    ## Not enough monotonicity constraints, completing with 0

    ## Warning in highs::hi_new_solver(model): WARNING: LP matrix packed vector
    ## contains 9 |value| in [8.88178e-16, 1.78572e-11] less than or equal to 1e-09:
    ## ignored

    ##    tau = 0.1 done

    ## Not enough monotonicity constraints, completing with 0

    ## Warning in highs::hi_new_solver(model): WARNING: LP matrix packed vector
    ## contains 9 |value| in [8.88178e-16, 1.78572e-11] less than or equal to 1e-09:
    ## ignored

    ##    tau = 0.5 done

    ## Not enough monotonicity constraints, completing with 0

    ## Warning in highs::hi_new_solver(model): WARNING: LP matrix packed vector
    ## contains 9 |value| in [8.88178e-16, 1.78572e-11] less than or equal to 1e-09:
    ## ignored

    ##    tau = 0.9 done
    ##
    ## > # Compute derivatives for analysis
    ## > compute_derivative <- function(y, x) {
    ## +   diff(y) / diff(x)
    ## + }
    ##
    ## > deriv_mixed <- compute_derivative(y_mixed, years_eval)
    ##
    ## > deriv_full <- compute_derivative(y_full_inc, years_eval)
    ##
    ## > # Find indices for analysis
    ## > idx_1945 <- which.min(abs(years_eval - 1945))
    ##
    ## > idx_1970 <- which.min(abs(years_eval - 1970))
    ##
    ## > cat("\n=== Derivative Analysis ===\n")
    ##
    ## === Derivative Analysis ===
    ##
    ## > cat(sprintf("Full increasing:       min derivative = %.4f (should be >= 0)\n",
    ## +             min(deriv_full)))
    ## Full increasing:       min derivative = -0.0000 (should be >= 0)
    ##
    ## > cat(sprintf("Mixed constraints:     min derivative = %.4f\n", min(deriv_mixed)))
    ## Mixed constraints:     min derivative = 0.0000
    ##
    ## > cat(sprintf("                       derivative between 1945-1970: range [%.4f, %.4f]\n",
    ## +             min(deriv_mixed[idx_1945:idx_1970]),
    ## +             max(deriv_mixed[idx_1945:idx_1970])))
    ##                        derivative between 1945-1970: range [0.0000, 0.0002]
    ##
    ## > # ============================================================
    ## > # Visualization
    ## > # ============================================================
    ## > # Ajouter un espace en haut de la page pour le titre
    ## > par(mfrow = c(2, 2), mar = c(4, 4, 4, 2), oma = c(0, 0, 2, 0))
    ##
    ## > # Plot 1: Unconstrained fit
    ## > plot(years, temperature, pch = 1,  col = "black",
    ## +      xlab = "Year", ylab = "Temperature Anomaly (C)",
    ## +      main = "1. Unconstrained Median (tau = 0.5)")

    ##
    ## > lines(years_eval, y_uncon, col = "red", lwd = 2)
    ##
    ## > abline(v = knots_years, col = "blue", lty = 2, lwd = 0.5)
    ##
    ## > abline(v = c(1945, 1970), col = "orange", lty = 3, lwd = 1.5)
    ##
    ## > abline(h = 0, col = "black", lty = 3, lwd = 0.5)
    ##
    ## > grid()
    ##
    ## > # Plot 2: Full monotonicity (too restrictive)
    ## > plot(years, temperature, pch = 1,  col = "black",
    ## +      xlab = "Year", ylab = "Temperature Anomaly (C)",
    ## +      main = "2. Full Increasing Constraint (too restrictive)")

    ##
    ## > lines(years_eval, y_full_inc, col = "blue", lwd = 2)
    ##
    ## > abline(v = knots_years, col = "blue", lty = 2, lwd = 0.5)
    ##
    ## > abline(v = c(1945, 1970), col = "orange", lty = 3, lwd = 1.5)
    ##
    ## > abline(h = 0, col = "black", lty = 3, lwd = 0.5)
    ##
    ## > text(1948, -0.1, "Forced increasing\nbut data shows cooling", col = "blue", cex = 0.7)
    ##
    ## > grid()
    ##
    ## > # Plot 3: Mixed constraints (decreasing 1945-1970, increasing elsewhere)
    ## > plot(years, temperature, pch = 1,  col = "black",
    ## +      xlab = "Year", ylab = "Temperature Anomaly (C)",
    ## +      main = "3. Mixed Constraints")

    ##
    ## > lines(years_eval, y_mixed, col = "darkgreen", lwd = 2)
    ##
    ## > abline(v = knots_years, col = "blue", lty = 2, lwd = 0.5)
    ##
    ## > abline(v = c(1945, 1970), col = "orange", lty = 2, lwd = 2)
    ##
    ## > abline(h = 0, col = "black", lty = 3, lwd = 0.5)
    ##
    ## > # Add shaded region for decreasing constraint
    ## > rect(1945, -0.8, 1970, 0.6, col = rgb(1, 0.5, 0, 0.1), border = NA)
    ##
    ## > text(1957, 0.5, "Decreasing constraint", col = "orange", cex = 0.8)
    ##
    ## > text(1900, 0.4, "Increasing", col = "darkgreen", cex = 0.7)
    ##
    ## > text(1985, 0.4, "Increasing", col = "darkgreen", cex = 0.7)
    ##
    ## > grid()
    ##
    ## > # Plot 4: Multiple quantiles with mixed constraints
    ## > colors <- c("orange", "darkgreen", "darkred")
    ##
    ## > plot(years, temperature, pch = 1,  col = "black",
    ## +      xlab = "Year", ylab = "Temperature Anomaly (C)",
    ## +      main = "4. Quantile Regression with Mixed Constraints")

    ##
    ## > for (i in seq_along(tau_multi)) {
    ## +   y_fit <- spline_eval(fits_multi[[i]], x_eval)
    ## +   lines(years_eval, y_fit, col = colors[i], lwd = 2,
    ## +         lty = ifelse(tau_multi[i] == 0.5, 1, 2))
    ## + }
    ##
    ## > abline(v = knots_years, col = "blue", lty = 2, lwd = 0.5)
    ##
    ## > abline(v = c(1945, 1970), col = "orange", lty = 2, lwd = 2)
    ##
    ## > abline(h = 0, col = "black", lty = 3, lwd = 0.5)
    ##
    ## > rect(1945, -0.8, 1970, 0.6, col = rgb(1, 0.5, 0, 0.1), border = NA)
    ##
    ## > legend("topleft", legend = c("tau = 0.1", "tau = 0.5", "tau = 0.9"),
    ## +        col = colors, lty = c(2, 1, 2), lwd = 2, cex = 0.8)
    ##
    ## > grid()
    ##
    ## > # Ajouter le titre global
    ## > mtext(paste("Global Temperature Anomaly Study 2 - Degree : ", degree),
    ## +       side = 3, line = 0.5, outer = TRUE, cex = 1.2, font = 2)
    ##
    ## > # ============================================================
    ## > # Conclusion
    ## > # ============================================================
    ## > cat("\n========================================\n")
    ##
    ## ========================================
    ##
    ## > cat("Conclusion\n")
    ## Conclusion
    ##
    ## > cat("========================================\n")
    ## ========================================
    ##
    ## > cat("The temperature data shows:\n")
    ## The temperature data shows:
    ##
    ## > cat("  - A general warming trend, especially after 1970\n")
    ##   - A general warming trend, especially after 1970
    ##
    ## > cat("  - A cooling period between 1945 and 1970\n")
    ##   - A cooling period between 1945 and 1970
    ##
    ## > cat("  - Full increasing constraint ignores this cooling period\n")
    ##   - Full increasing constraint ignores this cooling period
    ##
    ## > cat("  - Mixed constraints (decreasing 1945-1970, increasing elsewhere)\n")
    ##   - Mixed constraints (decreasing 1945-1970, increasing elsewhere)
    ##
    ## > cat("    better reflect the actual pattern\n")
    ##     better reflect the actual pattern
    ##
    ## > cat("  - Quantile regression reveals changing distribution over time\n")
    ##   - Quantile regression reveals changing distribution over time
    ##
    ## > cat("\nThis demonstrates that flexible constraint patterns\n")
    ##
    ## This demonstrates that flexible constraint patterns
    ##
    ## > cat("are more realistic than global monotonicity.\n")
    ## are more realistic than global monotonicity.
    ##
    ## > cat("\nDemo completed.\n")
    ##
    ## Demo completed.
    ##
    ## > par(oldpar)

    # Run specific demos:
    # demo("monotonicity")      - Increasing/decreasing constraints
    # demo("convexity")         - Convexity/concavity constraints
    # demo("der3")              - Third derivative constraints
    # demo("comprehensive")     - All constraint types
    # demo("temperature")       - Real-world example with temperature data

## Uniform Constraints

### Monotonicity Constraints

Monotonicity constraints force the fitted function to be non-decreasing (`monot = 1`) or non-increasing (`monot = -1`).

    # Generate data
    set.seed(42)
    n <- 100
    x <- seq(0, 1, length.out = n)
    y <- 2*x + 0.5*sin(4*pi*x) + 0.1*rnorm(n)
    knots <- quantile(x, probs = seq(0, 1, length.out = 10))

    # Unconstrained
    fit_uncon <- quantile_spline(x, y, knots, tau = 0.5)

    # Increasing constraint
    fit_inc <- quantile_spline(x, y, knots, tau = 0.5, monot = 1)

    # Decreasing constraint
    fit_dec <- quantile_spline(x, y, knots, tau = 0.5, monot = -1)

    # Compare
    x_eval <- seq(0, 1, length.out = 200)
    plot(x, y, pch = 16, cex = 0.5, col = "gray", main = "Monotonicity Constraints")
    lines(x_eval, spline_eval(fit_uncon, x_eval), col = "red", lwd = 2, lty = 2)
    lines(x_eval, spline_eval(fit_inc, x_eval), col = "blue", lwd = 2)
    lines(x_eval, spline_eval(fit_dec, x_eval), col = "green", lwd = 2)
    legend("topleft", legend = c("Unconstrained", "Increasing", "Decreasing"),
           col = c("red", "blue", "green"), lty = c(2, 1, 1), lwd = 2)

### Convexity Constraints

Convexity constraints force the second derivative to be non-negative (`convcons = 1`) or non-positive (`convcons = -1`).

    # Convex function
    y_conv <- (x - 0.5)^2 + 0.05*rnorm(n)

    fit_convex <- quantile_spline(x, y_conv, knots, tau = 0.5, convcons = 1)

    plot(x, y_conv, pch = 16, cex = 0.5, col = "gray", main = "Convexity Constraint")
    lines(x_eval, spline_eval(fit_convex, x_eval), col = "blue", lwd = 2)

### Third Derivative Constraints

Third derivative constraints control the curvature of the second derivative. The control of its sign allows enforcing smoothness of the curvature.

    # Cubic function with varying third derivative
    y_cubic <- x^3 - 1.5*x^2 + 0.5*x + 0.05*rnorm(n)

    fit_der3_pos <- quantile_spline(x, y_cubic, knots, tau = 0.5, der3cons = 1)

    plot(x, y_cubic, pch = 16, cex = 0.5, col = "gray", main = "Third Derivative Constraint")
    lines(x_eval, spline_eval(fit_der3_pos, x_eval), col = "blue", lwd = 2)

### Constraints by Degree

Degree | Monotonicity | Convexity | Third Derivative
--- | --- | --- | ---
1 (Linear) | +(linear) | X | X
2 (Quadratic) | + (Karlin) | +(constant) | X
3 (Cubic) | +(Karlin) | +(Karlin) | +(linear (constant))
4 (Quartic) | +(Karlin) | +(Karlin) | +(linear)

‘Karlin’ means in fact ‘quadratic’ inequality.

## Partial Constraints

Partial constraints apply monotonicity or convexity only on specific intervals, not the entire domain.

### Partial Monotonicity

    # Increasing only on first half
    monot_partial <- c(rep(1, 5), rep(0, 4))  # 5 intervals increasing, 4 unconstrained

    fit_partial <- quantile_spline(x, y, knots, tau = 0.5, monot = monot_partial)

    plot(x, y, pch = 16, cex = 0.5, col = "gray", main = "Partial Monotonicity")
    lines(x_eval, spline_eval(fit_partial, x_eval), col = "blue", lwd = 2)
    abline(v = knots[6], col = "red", lty = 2)  # Boundary of constrained region

### Partial Convexity

    # Convex on right half only
    conv_partial <- rep(0, length(knots))
    for (i in 1:length(knots)) {
      if (knots[i] > 0.5) conv_partial[i] <- 1
    }

    fit_conv_partial <- quantile_spline(x, y, knots, tau = 0.5, convcons = conv_partial)
    plot(x, y, pch = 16, cex = 0.5, col = "gray", main = "Partial Convexity")
    lines(x_eval, spline_eval(fit_conv_partial, x_eval), col = "blue", lwd = 2)
    abline(v = knots[6], col = "red", lty = 2)  # Boundary of constrained region

## Mixed Constraints

Multiple shape constraints can be applied simultaneously.

    # Increasing AND convex
    fit_both <- quantile_spline(x, y, knots, tau = 0.5,
                                monot = 1, convcons = 1)

    plot(x, y, pch = 16, cex = 0.5, col = "gray",
         main = "Mixed Constraints: Increasing + Convex")
    lines(x_eval, spline_eval(fit_both, x_eval), col = "blue", lwd = 2)

## Constraints with Multiple Knots

When using knot multiplicities, constraints need to be adapted to the reduced regularity.

    # Create knots with multiplicity at a point
    sn_mult <- c(0, 0, 0, 0, 0.3, 0.5, 0.5,0.5,0.5, 0.7, 1, 1, 1, 1)
    basis_mult <- Bspline_base(sn_mult, degree = 3)
    knots_mult <- basis_mult$knot

    # Fit with constraints
    fit_mult <- quantile_spline(x, y, knots_mult, tau = 0.5, monot = 1)
    plot(x, y, pch = 16, cex = 0.5, col = "gray",
         main = "Constraints with Multiple knots")
    lines(x_eval, spline_eval(fit_mult, x_eval), col = "blue", lwd = 2)

## Visualizing Constraint Effects

    par(mfrow = c(2, 2))

    # No constraints
    plot(x, y, pch = 16, cex = 0.4, col = "gray", main = "Unconstrained")
    lines(x_eval, spline_eval(fit_uncon, x_eval), col = "red", lwd = 2)

    # Monotonicity
    plot(x, y, pch = 16, cex = 0.4, col = "gray", main = "Monotonicity")
    lines(x_eval, spline_eval(fit_inc, x_eval), col = "blue", lwd = 2)

    # Convexity
    plot(x, y, pch = 16, cex = 0.4, col = "gray", main = "Convexity")
    lines(x_eval, spline_eval(fit_convex, x_eval), col = "blue", lwd = 2)

    # Mixed
    plot(x, y, pch = 16, cex = 0.4, col = "gray", main = "Mixed")
    lines(x_eval, spline_eval(fit_both, x_eval), col = "blue", lwd = 2)

    par(mfrow = c(1, 1))

## Karlin-Studden SOCP Formulation

The shape constraints are implemented using the Karlin-Studden characterization of non-negative polynomials of degree 2 or 3.

Depending on the degree of the spline (3 or 4), these constraints apply to the first or second derivative:

Spline Degree | Constraint Type | Polynomial Degree | Applied to
--- | --- | --- | ---
3 (Cubic) | Monotonicity | Quadratic (degree 2) | First derivative
3 (Cubic) | Convexity | Linear (degree 1) | Second derivative
4 (Quartic) | Monotonicity | Cubic (degree 3) | First derivative
4 (Quartic) | Convexity | Quadratic (degree 2) | Second derivative

Other constraints (third derivative, constraints at knots) are linear and are added as linear inequalities.

The complete set of constraints forms a Second-Order Cone Program (SOCP). This SOCP formulation ensures exact shape constraints, unlike approximations that only enforce constraints at knots or at a finite set of points.

The particular structure of the problem allows the use of the CVXR package, which handles this type of mixed quadratic/linear conic problems efficiently through its DCP (Disciplined Convex Programming) framework and supports multiple solvers (CLARABEL, OSQP, ECOS, SCS).

## Practical Considerations

### When to Use Shape Constraints

Constraint | When to Use
--- | ---
Monotonicity | Demand curves, growth curves, dose-response
Convexity | Cost functions, risk aversion, production functions
Third derivative | Smoothing curvature, spline regularity
Partial | Data with known local behavior

### Constraint Strength

Constraints can be too strong if: - The data clearly violates the assumed shape - Too few knots leading to over-smoothing - Multiple constraints conflicting

### Choosing the Degree

Degree | Smoothness | Flexibility
--- | --- | ---
1 | Piecewise linear | Low
2 | C¹ | Medium
3 | C² | High
4 | C³ | Very high

## Summary

The package provides: - Uniform constraints: Apply globally - Partial constraints: Apply on specific intervals - Mixed constraints: Multiple constraints simultaneously - SOCP formulation: Exact (not approximated) constraints - multiple knots: for lower regularity handling at knots (shocks)

For more examples, see the demos:

    demo(package = "BsplineQuantReg")

```

## Additional Sections You Could Add

Section | Content
--- | ---
Verification | How to check if constraints are satisfied
Comparison with cobs | Differences in constraint implementation
Real Examples | Temperature, economics, biology applications
Troubleshooting | Common issues and solutions