shape_constraints.knit

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