title: Force GAM model fit to be monotonic and go through a fixed point (x0, y0) with R mgcv
body: <p>I am trying to fit a GAM model to data under two constraints <em>simultatenously</em>: (1) the fit is monotonic (increasing), (2) the fit goes through a fixed point, say, <code>(x0,y0)</code>.</p>
  <p>So far, I managed to have these two constraints work <em>separately</em>:</p>
  <ul>
  <li><p>For (1), based on <a href="https://www.rdocumentation.org/packages/mgcv/versions/0.6-1/topics/pcls" rel="nofollow noreferrer">mgcv::pcls() documentation examples</a>, by using <code>mgcv::mono.con()</code> to get linear constraints sufficient for monotonicity, and estimate model coefs via <code>mgcv::pcls()</code>, using the constraints.</p>
  </li>
  <li><p>For (2), based on <a href="https://stat.ethz.ch/pipermail/r-help/2013-March/350253.html" rel="nofollow noreferrer">this post</a>, by setting the value of spline at knot location x0 to 0 + using offset term in the model formula.</p>
  </li>
  </ul>
  <p><strong>However, I struggle to combine these two constraints simultaneously.</strong> I <em>guess</em> a way to go is <code>mgcv::pcls()</code>, but I could work out neither (a) doing a similar trick of setting the value of spline at knot location x0 to 0 + using offset nor (b) setting equality constraint(s) (which I think could yield my (2) constraint setup).</p>
  <p>I also note that the approach for setting the value of spline at knot location x0 to 0 for my constrain condition (2) yields weirdly wiggly outcome (as compared to unconstrained GAM fit) -- as showed below.</p>
  <h3>Attempt so far: fit a smooth function to data under two constraints <em>separately</em></h3>
  <p>Simulate some data</p>
  <pre><code>library(mgcv)
  set.seed(1)
  x &lt;- sort(runif(100) * 4 - 1)
  f &lt;- exp(4*x)/(1+exp(4*x))
  y &lt;- f + rnorm(100) * 0.1
  dat &lt;- data.frame(x=x, y=y)
  </code></pre>
  <p>GAM unconstrained (for comparison)</p>
  <pre><code>k &lt;- 13
  fit0 &lt;- gam(y ~ s(x, k = k, bs = &quot;cr&quot;), data = dat)
  # predict from unconstrained GAM fit
  newdata &lt;- data.frame(x = seq(-1, 3, length.out = 1000))
  newdata$y_pred_fit0 &lt;- predict(fit0, newdata = newdata)
  </code></pre>
  <p><strong>GAM constrained: (1)</strong> the fit is monotonic (increasing)</p>
  <pre><code>k &lt;- 13
  
  # Show regular spline fit (and save fitted object)
  f.ug &lt;- gam(y~s(x,k=k,bs=&quot;cr&quot;))
  
  # explicitly construct smooth term's design matrix
  sm &lt;- smoothCon(s(x,k=k,bs=&quot;cr&quot;),dat,knots=NULL)[[1]]
  # find linear constraints sufficient for monotonicity of a cubic regression spline
  # it assumes &quot;cr&quot; is the basis and its knots are provided as input
  F &lt;- mono.con(sm$xp)
  
  G &lt;- list(
    X=sm$X,
    C=matrix(0,0,0),  # [0 x 0] matrix (no equality constraints)
    sp=f.ug$sp,       # smoothing parameter estimates (taken from unconstrained model)
    p=sm$xp,          # array of feasible initial parameter estimates
    y=y,           
    w= dat$y * 0 + 1  # weights for data    
  )
  G$Ain &lt;- F$A        # matrix for the inequality constraints
  G$bin &lt;- F$b        # vector for the inequality constraints
  G$S &lt;- sm$S         # list of penalty matrices; The first parameter it penalizes is given by off[i]+1
  G$off &lt;- 0          # Offset values locating the elements of M$S in the correct location within each penalty coefficient matrix.  (Zero offset implies starting in first location)
  
  p &lt;- pcls(G);       # fit spline (using smoothing parameter estimates from unconstrained fit)
  
  # predict 
  newdata$y_pred_fit2 &lt;- Predict.matrix(sm, data.frame(x = newdata$x)) %*% p
  # plot
  plot(y ~ x, data = dat)
  lines(y_pred_fit0 ~ x, data = newdata, col = 2, lwd = 2)
  lines(y_pred_fit2 ~ x, data = newdata, col = 4, lwd = 2)
  </code></pre>
  <p>Blue line: constrained; red line: unconstrained</p>
  <p><a href="https://i.stack.imgur.com/mC4aE.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/mC4aE.png" alt="enter image description here" /></a></p>
  <p><strong>GAM constrained: (2)</strong> fitted go through <code>(x0,y0)=(-1, -0.1)</code></p>
  <pre><code>k &lt;- 13
  
  ## Create a spline basis and penalty
  ## Make sure there is a knot at the constraint point (here: -1)
  knots &lt;- data.frame(x = seq(-1,3,length=k)) 
  # explicit construction of a smooth term in a GAM
  sm &lt;- smoothCon(s(x,k=k,bs=&quot;cr&quot;), dat, knots=knots)[[1]]
  
  ## 1st parameter is value of spline at knot location -1, set it to 0 by dropping
  knot_which &lt;- which(knots$x == -1)
  X &lt;- sm$X[, -knot_which]                  ## spline basis
  S &lt;- sm$S[[1]][-knot_which, -knot_which]  ## spline penalty
  off &lt;- dat$y * 0 + (-0.1)                   ## offset term to force curve through (x0, y0)
  
  ## fit spline constrained through (x0, y0)
  gam_1 &lt;- gam(y ~ X - 1 + offset(off), paraPen = list(X = list(S)))
  
  # predict (add offset of -0.1)
  newdata_tmp &lt;- Predict.matrix(sm, data.frame(x = newdata$x))
  newdata_tmp &lt;- newdata_tmp[, -knot_which]    
  newdata$y_pred_fit1 &lt;- (newdata_tmp %*% coef(gam_1))[, 1] + (-0.1)
  
  # plot
  plot(y ~ x, data = dat)
  lines(y_pred_fit0 ~ x, data = newdata, col = 2, lwd = 2)
  lines(y_pred_fit1 ~ x, data = newdata, col = 3, lwd = 2)
  # lines at cross of which the plot should go throught
  abline(v=-1, col = 3); abline(h=-0.1, col = 3)
  </code></pre>
  <p>Green line: constrained; red line: unconstrained</p>
  <p><a href="https://i.stack.imgur.com/F4ZGQ.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/F4ZGQ.png" alt="enter image description here" /></a></p>