Alexandru Ioan Cuza University of Iasi, Romania and
<EMAIL_ADDRESS>Ioan Cuza University of Iasi,
<EMAIL_ADDRESS>Andrei Arhire and Paul Diac [300]Theory of
computation Randomized local search [300]Theory of computation Graph
algorithms analysis The source code is available on Zenodo
(10.5281/zenodo.6646187) and GitHub
(https://github.com/AndreiiArhire/PACE2022).John Q. Open and Joan R. Access 2
42nd Conference on Very Important Topics (CVIT 2016) CVIT 2016 CVIT 2016
December 24–27, 2016 Little Whinging, United Kingdom 42 23

# PACE Solver Description: A Heuristic Directed Feedback Vertex Set Problem
Algorithm

Andrei Arhire  Paul Diac

###### Abstract

A feedback vertex set of a graph is a set of nodes with the property that
every cycle contains at least one vertex from the set i.e. the removal of all
vertices from a feedback vertex set leads to an acyclic graph. In this short
paper, we describe the algorithm for finding a minimum directed feedback
vertex set used by the _UAIC_ANDREIARHIRE_ solver, submitted to the heuristic
track of the 2022 PACE challenge.

###### keywords:

directed feedback vertex set, local search

###### category:

##  1 Preliminaries

Let $G=(V,E)$ be a directed graph.

  * •

$N^{-}_{G}(v)=\\{u\in V\mid\exists(u,v)\in E\wedge\nexists(v,u)\in E\\}$,
$d^{-}_{G}(v)=|N^{-}_{G}(v)|$;

  * •

$N^{+}_{G}(v)=\\{u\in V\mid\nexists(u,v)\in E\wedge\exists(v,u)\in E\\}$,
$d^{+}_{G}(v)=|N^{+}_{G}(v)|$;

  * •

$N^{\pm}_{G}(v)=\\{u\in V\mid\exists(u,v)\in E\wedge\exists(v,u)\in E\\}$,
$d^{\pm}_{G}(v)=|N^{\pm}_{G}(v)|$;

  * •

$G-v=(V\backslash\\{v\\},E\cap\\{\\{V\backslash\\{v\\}\\}\times\\{V\backslash\\{v\\}\\}\\})$;

  * •

$G\circ v$ =
$(V\backslash\\{v\\},E\cap\\{\\{\\{V\backslash\\{v\\}\\}\times\\{V\backslash\\{v\\}\\}\\}\cup\\{\\{N^{-}_{G}(v)\cup
N^{\pm}_{G}(v)\\}\times\\{N^{+}_{G}(v)\cup N^{\pm}_{G}(v)\\}\\}\\})$;

  * •

$G$ is a diclique if $\forall x,y\in V,x\neq y$ and $\\{(x,y),(y,x)\\}\subset
E$.

We will refer to $G-v$ operation as vertex $v$ removal (remove $v$ together
with all its adjacent edges) and to $G\circ v$ operation as vertex $v$ merger
(connect all its predecessors with all its successors and exclude $v$ from the
graph).

##  2 Solver Summary

The algorithm used by our solver has three main stages. In the first stage, it
finds an initial solution. The second stage consists of removing redundant
vertices contained in the initial solution. In the third stage, several local
searches are performed based on the best-known solution.

##  3 Reduction Rules

The problem can be viewed in the following way. We have to place each vertex
in a set A or a set B such that at the end, A is a minimum feedback vertex set
and B is the acyclic remainder. A group of vertices can be placed in A or B
based on some verification which can be done in a polynomial time.

The following two operations are considered to be reductions:

  * •

If a vertex can be part of set A without affecting the solution’s optimality,
we remove it and eventually introduce it in the solution set.

  * •

If a vertex can be part of set B without affecting the solution’s optimality,
then we merge it.

We make usage of 8 reduction rules described in [4], [2] and [3].

######  Reduction Rule .

If there exists a vertex $v\in V$ and an edge $(v,v)\in E$, remove $v$.

In this case, vertex $v$ contains a self-loop. It is erased and inserted into
the solution.

######  Reduction Rule .

If there exists a vertex $v\in V$, $(v,v)\not\in E$ with
$d_{G}^{-}(v)+d_{G}^{\pm}(v)\leq 1$ or $d_{G}^{+}(v)+d_{G}^{\pm}(v)\leq 1$,
merge $v$.

Vertex $v$ can have $d_{G}^{+}(v)=0$ or $d_{G}^{-}(v)=0$, then it can not be
part of any cycle. Otherwise, $d_{G}^{+}(v)=1$ or $d_{G}^{-}(v)=1$, connecting
all its predecessors with all its successors and excluding $v$ from the graph
(merge operation) will not affect the optimality. In both scenarios $v$ is not
considered in the solution.

######  Reduction Rule .

If there exists a vertex $v\in V$, $(v,v)\not\in E$,
$\min{(d_{G}^{-}(v),d_{G}^{+}(v))}=0$ and $N_{G}^{\pm}(v)$ forms a diclique,
remove $N^{\pm}_{G}(v)$ and merge $v$.

The critical observation here is that in a diclique, at most, one node will
not be part of the solution because any two vertices form a cycle. Thus, the
diclique can be viewed as a single node, so rule 2 can be applied further.

######  Reduction Rule .

If there exist vertices $u$, $v\in V$, $(u,v)\in E\wedge(v,u)\not\in E$ and
$u$ and $v$ are in different strongly connected components of the graph
$(V,E\backslash\\{E\cap\\{(x,y),\forall x\in V,y\in N^{\pm}_{G}(x)\\}\\})$,
erase $(u,v)$.

Edges between vertices in different strongly connected components are not part
of any cycle (otherwise, the vertices would be in the same component). Thus
all these edges can be deleted. Furthermore, in a diclique with two vertices,
at least one node will be removed, so edges that are part of a diclique can be
ignored when strongly connected components are computed.

######  Reduction Rule .

If there exist vertices $u$, $v\in V$, $(u,v)\in E\wedge(v,u)\not\in E$ and
$(N^{-}_{G}(u)\subset\\{N^{-}_{G}(v)\cup
N^{\pm}_{G}(v)\\})\vee(N^{+}_{G}(v)\subset\\{N^{+}_{G}(u)\cup
N^{\pm}_{G}(u)\\})$, erase $(u,v)$.

The idea with this rule is to delete a set of edges with the property that
there is no minimal cycle using them since only minimal cycles need to be
broken to compute the feedback vertex set.

The following three reduction rules are obtained based on the reduction rule 3
together with the idea that, at most, one node in a diclique is not part of
the solution.

######  Reduction Rule .

If there exists a node $v\in V$, $(v,v)\not\in E$ such that
$\\{N^{+}_{G}(v)\cup N^{\pm}_{G}(v)\\}$ or $\\{N^{-}_{G}(v)\cup
N^{\pm}_{G}(v)\\}$ forms a diclique, merge $v$.

######  Reduction Rule .

If there exists a node $v\in V$, $(v,v)\not\in E$ and $\\{N^{-}_{G}(v)\cup
N^{+}_{G}(v)\cup N^{\pm}_{G}(v)\\}$ can be split in two sets, $N^{\pm}_{G}(v)$
is included in one of them and each set forms a diclique, merge $v$.

######  Reduction Rule .

If there exists a node $v\in V$, $(v,v)\not\in E$, $d^{\pm}_{G}(v)=0$ and
$\\{N^{+}_{G}(v)\cup N^{-}_{G}(v)\\}$ can be split in at most three sets and
each set forms a diclique, merge $v$.

##  4 Stage one

In first part of the algorithm reductions rules are applied until the graph no
more supports any. If graph is not empty a promising vertex is selected to be
part of the solution and it is removed. These two operations are performed
until the graph becomes empty. A vertex $v$ is considered the best candidate
in the selection process if it maximizes among all the other vertices either
$(d^{+}_{G}(v)+d^{\pm}_{G}(v))\cdot(d^{-}_{G}(v)+d^{\pm}_{G}(v))$ or
$d^{\pm}_{G}(v)\cdot{\infty}+d^{-}_{G}(v)\cdot d^{+}_{G}(v)$.

##  5 Stage two

After many vertex selections, some of them may become redundant in the
solution, i.e., they are not part of any cycle so they can be excluded from
the feedback vertex set. To maximize the number of excluded vertices, we take
them in the reversed order of their insertion and introduce them in the
acyclic graph. At a fixed vertex, we check if the resulting graph is acyclic
or not using an optimized version of the Tarjan algorithm for computing
strongly connected components. If the graph is acyclic, the vertex will be
erased from the solution, and the graph will keep the changes. This step is
presented in [3] and [1].

##  6 Stage three

In the last part of the algorithm, some local searches are performed [3]. A
local search consists of running the same process as the one presented in
stage one but on a subgraph of G. Subgraphs are obtained by removing a
specific subset uniformly random from the best feedback vertex set found so
far together with a specific subset from the acyclic remainder. Additionally,
a slightly modified version of the algorithm from stage two is applied to the
new solution depending on the remaining time.

##  7 Complexity Analysis

Both selection criteria are simple but not always optimal, adding significant
benefits. Implying only the number of adjacent neighbors allows some
reductions to be performed efficiently using hash sets, keeping most of the
data in memory. After each update of the graph, reductions $1$ and $2$ can be
applied in $\mathcal{O}(1)$. Reduction $3$ is implemented in
$\mathcal{O}((|V|+|E|)\cdot\log{(|V|)})$ and reduction $4$ in
$\mathcal{O}(|V|+|E|)$. The rest are run only for vertices with degrees
bounded by a small constant, and we assume the complexity is
$\mathcal{O}(|V|+|E|)$. The first two reductions are run after each graph
change. The others are run after each loss of around $5\%-25\%$ of the number
of edges. For every iteration among $T$ ones, we choose to restore around 30%
of the vertices into the acyclic remainder and run the local search. As an
instance can be reduced in linear time to a graph with $|V|<|E|$, the final
complexity of the program is $\mathcal{O}(T\cdot(|V|+|E|\cdot\log{|E|}))$,
based on Master’s Theorem [5], where $T$ is the number of local searches. $T$
is a parameter bounded by the time limit.

## References

  * [1] Berend Hasselman.  An efficient method for detecting redundant feedback vertices.  CPB Discussion Paper 29, CPB Netherlands Bureau for Economic Policy Analysis, April 2004.  URL: https://ideas.repec.org/p/cpb/discus/29.html.
  * [2] Hen-Ming Lin; Jing-Yang Jou.  On computing the minimum feedback vertex set of a directed graph by contraction operations.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 19:295–307, 2000.  URL: http://doi.org/10.1109/43.833199, doi:10.1109/43.833199.
  * [3] Mile Lemaic.  Markov-Chain-Based Heuristics for the Feedback Vertex Set Problem for Digraphs.  PhD thesis, Universität zu Köln, 2008.
  * [4] Hanoch Levy; David W Low.  A contraction algorithm for finding small cycle cutsets.  Journal of Algorithms, 9:470–493, 1988.  URL: http://doi.org/10.1016/0196-6774%2888%2990013-2, doi:10.1016/0196-6774(88)90013-2.
  * [5] R.M. Verma.  A general method and a master theorem for divide-and-conquer recurrences with applications.  Journal of Algorithms, 16:67–79, 1994.  URL: http://doi.org/10.1006/jagm.1994.1004, doi:10.1006/jagm.1994.1004.


# On some extensions of shape-constrained generalized additive modelling in R

Natalya Pya Arnqvist
Department of Mathematics and Mathematical Statistics,
Umeå University, 901 87 Umeå, Sweden
<EMAIL_ADDRESS>
###### Abstract

Regression models that incorporate smooth functions of predictor variables to
explain the relationships with a response variable have gained widespread
usage and proved successful in various applications. By incorporating smooth
functions of predictor variables, these models can capture complex
relationships between the response and predictors while still allowing for
interpretation of the results. In situations where the relationships between a
response variable and predictors are explored, it is not uncommon to assume
that these relationships adhere to certain shape constraints. Examples of such
constraints include monotonicity and convexity. The `scam` package for R has
become a popular package to carry out the full fitting of exponential family
generalized additive modelling with shape restrictions on smooths. The paper
aims to extend the existing framework of shape-constrained generalized
additive models (SCAM) to accommodate smooth interactions of covariates,
linear functionals of shape-constrained smooths and incorporation of residual
autocorrelation. The methods described in this paper are implemented in the
recent version of the package `scam`, available on the Comprehensive R Archive
Network (CRAN).

Keywords: smoothing, shape constraints, interaction, smooth ANOVA, regression,
linear functionals of smooths

##  1 Introduction

Any application area that analyzes the relationship between a response and
multiple covariates could potentially benefit from using nonparametric and
semiparametric regression models. When analyzing the relationships between a
response and predictors, it might be natural to assume that some relationships
obey certain shape constraints, such as monotonicity and convexity. Such
problems are widespread in ecological and environmental studies. For example,
a county’s number of protest events is expected to increase with population
size. The relationship between tree height and altitude is known to decrease,
but tree height increases with tree age. The dose-effect curve in medicine,
the relationship between daily mortality and air pollution concentration, body
mass index and incidence of heart diseases are other examples where a shape
constraint is required. In such studies, unconstrained flexible nonparametric
modelling might be too flexible and give implausible or un-interpretable
results. Shape-constrained generalized additive models (SCAM), proposed in Pya
and Wood [2015] and implemented in an R package `scam` [Pya, 2024], provide a
general and powerful framework for semiparametric regression modelling with
shape constraints. SCAM is used in various application areas, including
ecological and environmental studies, education, public health, forestry,
genetics, medical and psychological research.

The additive models in the presence of shape constraints have been studied in
the larger literature on regression. Construction of shape-constrained spline-
based smoothers can be done by using linear inequality constraints on the
spline coefficients [Villalobos and Wahba, 1987, Ramsay, 1988, Kelly and Rice,
1990, Zhang, 2004, Wang and Xue, 2015, Meyer, 2018]. However, the linear
inequality constraints bring methodological difficulties in estimating the
degree of spline smoothness. Several other papers have studied additive shape-
constrained regression based on B-splines, including He and Shi [1998],
Bollaerts et al. [2006], Rousson [2008], Wang and Meyer [2011], Wang and Zhou
[2021]. Groeneboom et al. [2001, 2008], Guntuboyina and Sen [2015] has
considered the problem of univariate convex regression. Alternative approaches
are based on non-spline ideas. Mammen [1991], Hall and Huang [2001], Du et al.
[2013] focused on shape constrained kernel regression techniques,
whereasAntoniadis et al. [2007] proposed a penalized wavelet regression. There
was also work in the direction of additive models without smoothing that has
been covered by Mammen and Yu [2007], Cheng [2009], Fang and Meinshausen
[2012], Cheng et al. [2012], Yu [2014], Chen and Samworth [2016]. The problem
of monotonic additive regression with Bayesian ideas was considered by Holmes
and Heard [2003], Lang and Brezger [2004], Dunson and Neelon [2003], Dunson
[2005], and with boosting approach by Tutz and Leitenstorfer [2007]. Brezger
and Steiner [2008] developed a Bayesian approach for GAM using penalized
B-splines and incorporated monotonicity assumption on univariate smooth by
introducing specific prior distributions. Despite such diverse literature on
shape-constrained regression, existing procedures either lack the ability to
efficiently estimate the smoothing parameters in multiple model components
setting or result in non-smooth functions.

The models with multi-dimensional smooths under shape constraints have been
much less studied. Bollaerts et al. [2006] introduced bivariate shape-
constrained P-splines in the least square setting with two covariates only.
Lin and Dunson [2014] proposed a Bayesian regression method for estimating
monotone multi-dimensional functions, building on Gaussian process projection.
Dette and Scheder [2006] developed monotone smooth regression for several
covariates based on kernel regression. A slightly more general case was
considered by Du et al. [2013], who extended the idea of Hall and Huang [2001]
of monotone kernel regression smoothing to a multivariate setting, covering
monotonicity and concavity. Except for the SCAM, the existing approaches work
only with normal responses and a non-additive model structure.

The main objective of this paper is to provide an extension of the existing
framework for generalized additive modelling with a mixture of unconstrained
terms and various shape-restricted terms to accommodate smooth interaction of
covariates, varying coefficient terms, linear functionals with or without
shape constraints as model components, and data with short-term temporal or
spatial autocorrelation. This paper contributes to the following novel
elements: i) bivariate interaction smooths with increasing and decreasing
constraints are proposed, which allows for functional ANOVA decomposition
(smooth ANOVA) within generalized additive modelling with shape constraints;
ii) extension of SCAM with AR1 model on the residuals is introduced, where
residual auto-correlation with an AR1 correlation structure can be dealt with
in Gaussian models with identity links; iii) extension of SCAM to allow linear
functionals of shape-constrained smooths as model components, is proposed,
which is known as scalar-on-function regression in the field of functional
data analysis.

The remainder of the paper is structured as follows. The general modelling
framework is introduced in Section 2. Section 3 describes the SCAM extensions
and illustrates their use in examples. Conclusions and some discussions are
given in Section 4.

##  2 The general SCAM framework

The class of models considered can be written as

$g(\mu_{i})={\bf
A}_{i}{\bm{\theta}}+\sum_{j}f_{j}(z_{ji})+\sum_{k}m_{k}(x_{ki}),\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ Y_{i}\sim{\rm EF}(\mu_{i},\phi),$
(1)

where $g$ is a specified link function, $Y_{i}$ is a random real-valued
response, for $i=1,\ldots,n,$ with the mean $\mu_{i}$ with an exponential
family distribution (scale parameter $\phi$), ${\bf A}_{i}$ is the $i$th row
of a known parametric model matrix with an associated vector of unknown
coefficients $\bm{\theta},$ $f_{j}$ is an unknown smooth function of covariate
$z_{j}$, and $m_{k}$ is an unknown shape constrained smooth function of
covariate $x_{k}$. The covariates $z_{j}$ and $x_{k}$ may be vector-valued.

What makes models in (1) different from the standard generalized additive
models (GAM) are the shape constraints on $m_{k}.$ These models gained
popularity through Pya and Wood [2015], which provided computational efficient
estimation methods. The work recognizes that the full practical benefits from
allowing flexible dependence on covariates could only be realized if the
degrees of smoothing of the $f_{j}$ and $m_{k}$ are estimated as part of model
fitting. Pya and Wood [2015] proposed shape-constrained P-splines (SCOP-
splines) based on a nonlinear reparametrization of P-splines [Eilers and Marx,
1996] with discrete penalties, which allows a variety of shape constraints for
univariate and multivariate smooths. Specifically, for a univariate smooth

$m(x)=\sum_{j=1}^{q}\gamma_{j}b_{j}(x),$

where $q$ is the number of basis functions, $b_{j}$ are B-spline basis
functions of at least second order, $\gamma_{j}$ are unknown coefficients.
Using the first order derivative feature of the B-splines, shape constraints
are imposed by re-parameterizing the spline coefficients
$\bm{\gamma}=\bm{\Sigma}\tilde{\bm{\beta}},$ where
$\bm{\gamma}=(\gamma_{1},\ldots,\gamma_{q})^{T},$ $\bm{\Sigma}$ is a matrix of
0’s and 1’s, $\bm{\beta}=(\beta_{1},\ldots,\beta_{q})^{T},$ and
$\tilde{\bm{\beta}}$ is a vector of $\bm{\beta}$ elements, where depending on
the shape constraints, some of the elements are exponentiated $\beta_{j}$. For
example, in case of a monotonically increasing smooth,
$\tilde{\bm{\beta}}=\left[\beta_{1},\exp(\beta_{2}),\ldots,\exp(\beta_{q})\right]^{T},$
while $\Sigma_{ij}=0$ if $i<j$ and $\Sigma_{ij}=1$ if $i\geq j.$ The $n$
vector of the evaluated function values, $m(x),$ can now be expressed as ${\bf
m}={\bf X}\bm{\Sigma}\tilde{\bm{\beta}},$ where ${\bf X}$ is a matrix with
$X_{ij}=b_{j}(x_{i}).$ In a smoothing context, with each $m(x)$ is associated
a smoothing penalty, which is quadratic in the $\beta_{j}$ and controls the
‘wiggliness’ of $m.$ In a matrix form, the penalty can be written as
$\bm{\beta}^{T}{\bf S}\bm{\beta},$ where $\bf S$ is a matrix of known
coefficients. Multivariate shape-constrained smooths with constraints assumed
on either all or a selection of the covariates are built up using the concept
of tensor product splines.

To represent (1) for computation, all unconstrained $f_{j}$ are represented
via basis expansions of modest rank, with chosen penalties and identifiability
constraints. This allows to specify $\sum_{j}f_{j}(z_{ji})={\bf
F}_{i}\bm{\gamma},$ where $\bf F$ is a model matrix, $\bm{\gamma}$ is a vector
of unknown coefficients. The penalties on the $f_{j}$ are quadratic in
$\bm{\gamma}.$ By absorbing the corresponding identifiability constraints, the
model matrices for all the $m_{k}$ can be combined so that
$\sum_{k}m_{k}(x_{ki})={\bf M}_{i}\tilde{\bm{\beta}},$ where $\bf M$ is a
model matrix, $\tilde{\bm{\beta}}$ is a vector containing a mixture of
$\beta_{j}$ and exponentiated $\beta_{j}.$ So (1) becomes $g(\mu_{i})={\bf
A}_{i}{\bm{\theta}}+{\bf F}_{i}\bm{\gamma}+{\bf M}_{i}\tilde{\bm{\beta}},$ and
by combining the matrices column-wise into one model matrix ${\bf
X}=\left[{\bf A}:{\bf F}:{\bf M}\right],$

$g(\mu_{i})={\bf X}_{i}\tilde{\bm{\beta}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ Y_{i}\sim{\rm EF}(\mu_{i},\phi),$
(2)

where $\tilde{\bm{\beta}}$ now contains $\bm{\theta},$ $\bm{\gamma}$ and
original $\tilde{\bm{\beta}}.$ Let $l(\bm{\beta})$ be the model log likelihood
and ${\bf S}^{k}$ is the penalty matrix of known coefficients that correspond
to the $k$th smoothing penalty. Generally, ${\bf S}^{k}$ has only a small
block of non-zero elements. The estimated model coefficients are then

$\hat{\bm{\beta}}=\underset{\bm{\beta}}{\textrm{argmax}}\left\\{l(\bm{\beta})-\frac{1}{2}\sum_{k}^{M}\lambda_{k}\bm{\beta}^{T}{\bf
S}^{k}\bm{\beta}\right\\},$ (3)

given $M$ values for the smoothing parameter $\lambda_{k},$ controlling the
level of penalization. Smoothing parameters are estimated to minimize a
prediction error criterion such as GCV or AIC in the outer optimization
scheme. Since the imposed extra non-linearity of SCOP-splines coefficients
bans re-using or modifying the methods developed for unconstrained GAMs,
different computational schemes for estimating the model coefficients,
$\bm{\beta},$ and smoothing parameters, $\lambda_{k},$ were developed, as well
as simulation-free approximate Bayesian confidence intervals for the model
terms Pya and Wood [2015]. Over the last few years, various SCOP-splines have
been introduced (upon scam-users’ requests) and implemented in the `scam`
package. The univariate single penalty built-in shape-constrained smooth
classes include such constraints as monotonicity, concavity, a mixture of
monotonicity and concavity, monotonicity plus the restriction on passing
through zero at the right-end point (left-end point) of the covariate range,
and positivity constraint. Seventeen bivariate SCOP-splines under monotonicity
and/or concavity constraints are suggested, where monotonicity (concavity) may
be assumed on only one of the covariates (single monotonicity) or both of them
(double monotonicity) are available as the smooth terms of the SCAM. Double or
single monotonicity (concavity) is achieved by the corresponding re-
parametrization of the bivariate basis coefficients to satisfy the sufficient
conditions formulated in terms of the first order differences of the
coefficients. Double penalties for the shape-constrained tensor product
smooths are obtained from the penalties of the marginal smooths.

##  3 Extended SCAM

To accommodate the extended SCAMs proposed in this paper, the class of models
in (1) can now be written as

$g(\mu_{i})={\bf
A}_{i}{\bm{\theta}}+\sum_{j}L_{ij}f_{j}+\sum_{k}L_{ik}m_{k},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ Y_{i}\sim{\rm EF}(\mu_{i},\phi),$
(4)

where $L_{ij}$ and $L_{ik}$ are bounded linear functionals, as in (1) $f_{j}$
and $m_{k}$ are unconstrained and shape-constrained smooth functions of some
known covariates correspondingly. In the standard case $L_{ik}$ (and $L_{ij}$)
is simply a functional of evaluation, $L_{ik}m_{k}=m_{k}(x_{i}),$ however for
varying coefficient models $L_{ik}m_{k}=m_{k}(x_{i})z_{i},$ for scalar on
function regression $L_{ik}m_{k}=\int s_{i}(x)m_{k}(x)dx.$ Also, a linear
random effects term ${\bf U}_{i}{\bf b},$ where ${\bf b}\leavevmode\nobreak\
\sim\leavevmode\nobreak\ N({\bf 0},{\bf I}\sigma_{u}^{2})$ and $\bf U$ is a
model matrix that can be added to the model structure of (4), thus extending
the SCAM to shape constrained generalized additive mixed models in the same
way as it is done for unconstrained GAM Wood [2017]. The following subsections
describe the smooth-ANOVA with shape-constrained smooths, shape-constrained
additive modelling with simple autocorrelation on the residuals, and linear
functional extensions for SCAM.

###  3.1 SCAM with smooth interactions

Considering models structured with main effects and interactions can be
beneficial in certain scenarios, for instance:

$f_{1}(x)+f_{2}(z)+m_{3}(x,z).$

Here, $f_{1}$ and $f_{2}$ represent smooth ‘main effects’. $m_{3}$ is a smooth
‘interaction’ subject to, for example, an increasing constraint with respect
to $x$ and the main effects excluded ($f_{1}$ can also be subject to an
increasing constraint). The construction of such a functional ANOVA
decomposition of shape-constrained smooths relies on constructing a tensor
product basis. This recognizes the equivalence between tensor product basis
construction and the method used for interactions in linear models. Regression
modelling with interactions within generalized additive models, where the will
is to separate the unconstrained smooth interaction effect of two or more
covariates from the main smooth effects (ANOVA decompositions of smooths), was
introduced in Wood [2017]. Following that work, the marginal smooths of a
tensor product are subject to identifiability constraints before constructing
the tensor product basis. This leads to interaction smooths that do not
encompass the corresponding main effects.

Two tensor product bivariate interaction smooths under shape constraints have
been developed and implemented in the R package `scam`, with increasing and
decreasing constraints in the first covariates. Both smooth interactions apply
the corresponding SCOP identifiability constraints to the first marginal and
centering constraints to the second unconstrained marginal. The SCOP
identifiability and centering constraints remove the constant functions from
the bases of the marginals. By doing so, constant functions are eliminated
from the marginal bases, preventing the interaction term from containing the
main effects of $f_{1}$ and $f_{2}$. This would happen if the marginal bases
were multiplied by the constant functions in the other marginal bases. The
removal of constraints does not impact penalty terms, given that constant
functions lie within the null space of the penalty. In the `scam` package, two
shape-constrained smooth interaction terms with increasing and decreasing
constraints are implemented.

As an illustration, consider the `wesdr` data set from the `gss` package of
Chong Gu. This data set originates from an epidemiological study of diabetic
retinopathy. It contains 669 observations across four variables: duration of
diabetes (`dur`, in years), percent of glycosylated hemoglobin in the blood
(`gly`), body mass index (`bmi`), and binary indicator of diabetic retinopathy
progression (`ret`). To explore and understand the data, we can start by
modelling an unconstrained smooth additive logistic regression model as
suggested in Wood [2020]. We might expect interactions between covariates, so
we include smooth interaction terms. The model structure becomes:

$\begin{split}\mathrm{logit}(p_{i})=&f_{1}(\mathtt{dur}_{i})+f_{2}(\mathtt{gly}_{i})+f_{3}(\mathtt{bmi}_{i})+\\\
&f_{4}(\mathtt{dur}_{i},\mathtt{gly}_{i})+f_{5}(\mathtt{dur}_{i},\mathtt{bmi}_{i})+f_{6}(\mathtt{gly}_{i},\mathtt{bmi}_{i}),\end{split}$

where $\mathtt{ret}_{i}\leavevmode\nobreak\ \sim\leavevmode\nobreak\
\mathrm{Binom}(1,p_{i}),$ $p_{i}$ is probability of retinopathy,
$f_{1},f_{2},f_{3}$ are smooth main effects, and $f_{4},f_{5},f_{6}$ are
smooth interaction terms. The estimated smooths of the unconstrained model
shown in figure 7 of Wood [2020] indicate a non-zero interaction effect
$f_{6}.$ However, we might expect that if there is such an interaction effect
of $\mathtt{gly}_{i}$ and $\mathtt{bmi}$, the risk of retinopathy will not
decrease with increasing body mass index. The unconstrained smooth ANOVA model
might be too flexible to catch such an effect. Increasing constraints on the
interaction smooths along the $\mathtt{bmi}$ covariate could result in more
plausible effects. So, a shape-constrained smooth ANOVA model can be
considered here, and it is easily fitted

    m <- scam(ret ~ s(dur,k=k)+ s(gly,k=k)+ s(bmi,k=k)+ ti(dur,gly,k=k)
             + s(bmi,dur,bs="tismi") + s(bmi,gly,bs="tismi"),
             data=wesdr, family=binomial())

The estimated smooths are shown in figure 1 and reveal sounder results for
$f_{6}$ with the non-decreasing risk effect along the body mass index.
Introducing shape-constrained interaction smooths leads to only a small
increase in AIC.

Figure 1: Estimated smooth terms of the shape-constrained smooth ANOVA
logistic regression model for diabetic retinopathy example.

###  3.2 SCAM with AR1 model on the residuals

Many environmental and business data sets show short-term temporal, spatial or
spatial-temporal autocorrelation. It appears that the fitting problem for non-
independent data with a general multivariate normal distribution with an
unknown mean and a covariance matrix known up to a constant of proportionality
can be transformed to match the fitting problem for independent data with a
constant variance exactly. An AR(1) error model can be added to SCAM to deal
with the autocorrelation used for the residuals of Gaussian-identity link
models in the same way suggested for unconstrained models Wood et al. [2015].
The model can be expressed as $Y_{i}={\bf
X}_{i}\tilde{\bm{\beta}}+\epsilon_{i},$ where
$\epsilon_{i}=\rho\epsilon_{i-1}+e_{i}$ and $e_{i}$ are independent
$N(0,\sigma^{2}).$ The covariance matrix for this model,
$\mathrm{cov}(\bm{\epsilon})=\sigma^{2}{\bf V}$ is

${\bf V}=\left(\begin{matrix}1&\rho&\rho^{2}&\ldots&\rho^{n-1}\\\
\rho&1&\rho&\ldots&\rho^{n-2}\\\ \ldots&\ldots&\ldots&\ldots&\ldots\\\
\rho^{n-1}&\rho^{n-2}&\rho^{n-3}&\ldots&1\end{matrix}\right),$

while the Cholesky factor of the inverse, ${\bf V}^{-1}={\bf C}^{T}{\bf C}$ is

${\bf C}=\frac{1}{\sqrt{1-\rho^{2}}}\left(\begin{matrix}1&-\rho&0&\ldots&0\\\
0&1&-\rho&\ldots&0\\\ \ldots&\ldots&\ldots&\ldots&\ldots\\\
0&0&0&\ldots&\sqrt{1-\rho^{2}}\end{matrix}\right),$

which is banded, having only two non-zero diagonals. So, if we set $\tilde{\bf
Y}={\bf CY}$ and $\tilde{\bf X}=\bf CX,$ we have $\tilde{\bf Y}=\tilde{\bf
X}\tilde{\bm{\beta}}+\tilde{\bm{\epsilon}},$ where $\tilde{\bm{\epsilon}}={\bf
C}\bm{\epsilon}$ are now independent $N(0,\sigma^{2})$. The resulting model is
in the form of (2), so the estimation methods and inference about
$\tilde{\bm{\beta}}$ developed for SCAM with independent data can be used
here. The banded structure of the Cholesky factor $\bf C$ makes the
transformation by $\bf C$ computationally cheap. Formation of $\tilde{\bf Y}$
and $\tilde{\bf X}$ involves simple weighted differencing operations on the
elements of $\bf Y$ and $\bf X$ with $O(n)$ and $O(np)$ operations
correspondingly. When $\rho$ is unknown, AIC or GCV can be used to estimate it
in a simple one-dimensional direct search. In this case, $\bf C$ should be
obtained for every trial $\rho.$ However, as $\bf C$ can be obtained without
forming ${\bf V}^{-1}$ and it has the banded structure, the suggested AR(1)
error model is computational feasible. If $\rho$ is supplied (as known), then
$\bf C$, $\tilde{\bf Y}$ and $\tilde{\bf X}$ are formed only once.

Consider `sitka` data set from an R package `SemiPar` containing log-size
measurements for 79 Sitka spruce trees grown in normal (control) or ozone-
enriched conditions. Figure 2 shows the growth trajectories of these trees. An
unconstrained additive mixed modelling for this data (without residual
autocorrelation) was illustrated in Wood [2017]. The model estimated smooth
mean growth curve reveals a dubious decreasing effect for times in days
between 400 and 500 (top left plot in figure 3), so it is reasonable to impose
an increasing constraint on the growth smooth term. In addition, the model
residuals contain some autocorrelation (top right plot in figure 3). So, an
AR1 model for the residuals can be included to reduce the effects of the
autocorrelation. A simple shape-constrained additive mixed model with an AR1
error model would be

$\mathtt{log.size}_{i}=m(\mathtt{day}_{i})+\mathtt{ozone}_{i}+b_{\mathtt{id}_{i}}+\epsilon_{i},\leavevmode\nobreak\
\leavevmode\nobreak\ b_{\mathtt{id}_{i}}\leavevmode\nobreak\
\sim\leavevmode\nobreak\ N(0,\sigma^{2}_{b}),$

where $\epsilon_{i}$ is an AR1 error model within each tree,
$m(\mathtt{day}_{i})$ is an increasing smooth term, $\mathtt{ozone}_{i}$ is a
binary factor of ozone treatment, $\mathtt{id}_{i}$ is the tree identification
number to which the $i^{\mathrm{th}}$ observation belongs. The following code
would estimate such a model:

    s <- scam(log.size ~ s(days,bs="mpi",k=15)+ozone+s(id.num,bs="re"),
             data=sitka, AR1.rho=rho, AR.start=start.event)

Here, `AR.start` argument is a logical vector indicating where the breaks
between autocorrelated sequences of residuals are needed, which is `TRUE` at
the first observation of each tree and `FALSE` otherwise. `AR1.rho` is the AR1
correlation parameter, which can be chosen by trying a sequence of values and
selecting the value that minimizes the AIC score. To get a likely range of
`rho` values, the autocorrelation of the residuals of the fitted un-correlated
model can be used. The lower right plot of figure 3 shows the estimated
autocorrelation function of the standardized residuals of the fitted model
returned in `s$std.rsd`. The ideal AR1 model would result in standardized
residuals that appear approximately uncorrelated. Although not perfect, the
model is largely improved compared to the uncorrelated model (top right plot),
resulting in a much worse ACF. The fitted model `s` showed a significant
negative effect of the enriched ozone, and figure 3 demonstrates the estimated
smooth mean Sitka spruce tree growth curve (lower left) with 95% credible
interval and the normal QQ-plot for the predicted random effect from the model
(lower middle).

Figure 2: Growth paths for 79 Sitka spruce trees: (left) paths for 25 trees
grown in control environments; (right) paths for 54 trees grown in ozone-
enriched environments. Figure 3: The top row is for the unconstrained
uncorrelated model while the lower row is for the scam with AR1 error model,
s. The first column plots the estimated smooth mean growth curve of Sitka
spruce trees with 95% credible interval. The next column plots the normal QQ-
plot for the predicted random effects for the models. The right column shows
the ACF plots of the residuals for the un-correlated model (top) and of the
standardized residuals for the model s (lower).

###  3.3 Linear functionals of shape constrained smooths

Since the model setup of SCAM is based on the model setup of a standard GAM of
the R package `mgcv` Wood [2023], the methodology developed for SCAM allows
for including any unconstrained smooth terms of one or more variables
available in the `mgcv` into its model structure. So, the linear functionals
of smooth functions introduced in GAM can be extended to shape-constrained
versions. One such example is a scalar on function regression, as known in the
field of functional data analysis, which describes the relationship between a
scalar response and a set of functional covariates. Suppose that, for subject
$i=1,\cdots,n$, we observe data of the form
$\left\\{y_{i},(z_{ij},t_{ij})\right\\},\leavevmode\nobreak\ j=1,\cdots,J,$
where $y_{i}$ is a scalar response and $z_{ij}$ are discrete realizations of a
real-valued, square-integrable smooth curve $z_{i}()$ at observation points
$t_{ij}$. A simple functional linear model can be written as

$Y_{i}=\alpha+\int m(t)z_{i}(t)dt+\epsilon_{i},\leavevmode\nobreak\
\leavevmode\nobreak\ \epsilon_{i}\leavevmode\nobreak\ \sim\leavevmode\nobreak\
N(0,\sigma^{2}),$

where $m(t)$ is a smooth function to estimate subject to some shape
constraints. Using the SCOP-spline basis expansion
$m(t)=\sum_{j=1}^{q}\gamma_{j}b_{j}(t),$ the model matrix elements of the
linear functional term become $X_{ij}=\int b_{j}(t)z_{i}(t)dt.$ Given the
linearity of the model in $m(t),$ and by replacing the integral with a
discrete sum approximation (numerical quadrature), the model estimation
follows the exact methodology developed for SCAMs. The model setup of the
`mgcv` is used to specify the model terms, where a summation convention
(summation across rows) is applied when a metric `by` variable and a smooth
covariate are expressed as matrices.

To fit the functional linear model with, for example, decreasing constraint,
the `scam` call would be

    b <- scam(y ~ s(X,by=Z,bs="mpdBy"))

The linear functional term is specified via `s()` smooth with a `by` variable,
where `Z` and `X` are matrices of the same dimensions, `Z` is a matrix of
discrete observations of the smooth functional predictors, `X` is a matrix of
the observation points, $t_{ij}.$ The type of shape constraints imposed on the
smooth coefficient function, $m(t)$, is specified in the `bs` argument.
`bs="mpdBy"` in the example above indicates that $m(t)$ is subject to
decreasing constraint. Eight different SCOP-splines with monotonicity,
concavity and a mixture of monotonicity and concavity constraints are
implemented in the `scam` package to be used with a numeric ‘by’ variable that
takes more than one value. The construction of these smooth terms is similar
to that of the ordinary SCOP-splines, with the only difference being that they
are built without applying identifiability constraints, since smooth terms
with a numeric ‘by’ variable taking more than one value are fully
identifiable. Figure 4 shows a simulated example of the linear functionals of
the smooth function under decreasing constraint. From the two lower right
plots, it is clear that the shape-constrained model performs better than the
unconstrained version.

Figure 4: Simulated example of the scalar on function regression with
decreasing constraint. The first six plots (from top left to lower right) show
six examples of the generated functional predictors. In total, 200 functions
were simulated and stored in the matrix Z, each at 100 evaluation points.
Third lower plot: estimated constrained coefficient function (solid line) and
true coefficient function (dashed line). Lower right plot: estimated
coefficient function for the unconstrained scalar on function regression
(solid line) and true coefficient function (dashed line).

Note that function on scalar regression is also possible with SCAM by
expanding the functional coefficients with unconstrained or SCOP-spline basis
expansion.

$Y_{i}(t)=f_{0}(t)+\sum\limits_{j}x_{ij}f_{j}(t)+\sum\limits_{k}x_{ik}m_{k}(t)+\epsilon_{i}(t),$

where $Y_{i}$ are functional responses, $x_{i}$ are scalar predictors, $f_{j}$
and $m_{k}$ are functional coefficients to estimate correspondingly without or
with shape constraints. Errors, $\epsilon_{i}(t)$, correlated within a subject
can be modelled as AR1 error process with `scam`.

###  3.4 Smoothing parameter estimation

The current fitting approach for SCAM is reasonably efficient and robust when
applied to data sets containing a few tens of thousands of observations. If
$n$ and $p$ denote numbers of rows and columns of the model matrix
correspondingly, then the approach has the $O(np^{2})$ computational cost of
standard penalized regression spline based GAM estimation and typically
involves 2-4 times as many $O(np^{2})$ steps due to the additional non-
linearities required for dealing with the shape constrained terms. As an
alternative to the conventional quasi-Newton approach for smoothing parameter
$\bm{\lambda}$ optimization, an extended Fellner-Schall method (EFS) reported
in Wood and Fasiolo [2017] was introduced. The smoothing parameters are
estimated in the outer optimization scheme using a full Newton method for
coefficient estimation. The extended Fellner-Schall method generalizes the
technique of Fellner [1986] and Schall [1991] and offers a significant
simplification over existing methods by re-using quantities anyway required
for model coefficients estimation when iteratively updating $\bm{\lambda}$.
When modelling with SCAM, rather than maximizing the log restricted marginal
likelihood as it is proposed in [Wood and Fasiolo, 2017], similar convenient
updates are obtained for minimizing a prediction error GCV criterion,
$\mathcal{V}_{g}=nD(\hat{\bm{\beta}})/(n-\gamma\tau)^{2},$ when the scale
parameter of the exponential family distribution is unknown and for the UBRE,
$\mathcal{V}_{a}=D(\hat{\bm{\beta}})+2\phi\tau\gamma,$ for known scale
parameter. The smoothing parameter update for the GCV minimization is

$\rho_{j}^{*}=\rho_{j}+\log\left\\{-\frac{2\gamma
D}{n-\tau}\cdot\frac{\partial\tau}{\partial\rho_{j}}\left[\left(\frac{\partial
D}{\partial\bm{\beta}}\right)^{T}\cdot\frac{d\hat{\bm{\beta}}}{d\rho_{j}}\right]^{-1}\right\\},$

where $\rho_{j}=\log(\lambda_{j})$ to ensure that $\lambda_{j}^{*}$ is
positive, $D$ is the model deviance, $\tau$ is the effective degrees of
freedom, $\gamma$ is a parameter that usually has the value $1,$ but can be
increased to get smoother models. The $\rho$ update when minimizing the UBRE
score is

$\rho_{j}^{*}=\rho_{j}+\log\left\\{-2\phi\gamma\cdot\frac{\partial\tau}{\partial\rho_{j}}\left[\left(\frac{\partial
D}{\partial\bm{\beta}}\right)^{T}\cdot\frac{d\hat{\bm{\beta}}}{d\rho_{j}}\right]^{-1}\right\\}.$

The EFS approach eliminates the need for implicit differentiation in quasi-
Newton optimization and leads to a more robust smoothing parameter scheme for
SCAM. It is implemented in the `scam` package and is called for by
`optimizer="efs"` argument of the `scam()` function.

##  4 Discussion

To date, all developed methods within shape-constrained regression modelling
are limited to either modelling with only independent normal response
variables or models with standard additive model structures of smooth
functions. The main contribution of this work has been to provide practical
extensions to the SCAM framework for fitting models that include shape-
constrained smooth interactions and linear functionals of the shape-
constrained smooths as model components and can deal with residual
autocorrelation for the Gaussian-identity link models. The proposed methods
seamlessly integrate into the existing general framework of shape-constrained
generalized additive modelling, SCAM. So, anyone familiar with generalized
additive modelling and shape-constrained additive modelling through R packages
such `mgcv` and `scam` can readily apply these extensions. The methods are
implemented in the R package `scam`, which is under continual development and
maintenance.

The implementation cost of the suggested EFS method for estimating smoothing
parameters is significantly reduced, thus allowing for easier application and
opening doors for a wider range of models with shape-constrained terms. The
costs can be further reduced by applying a quasi-Newton method for model
coefficients estimation in which the Hessian of the log-likelihood will be
replaced in the update by a quasi-Newton approximation, thus only requiring
the gradient of the log-likelihood. This approach is added as an alternative
fitting method (working with the EFS for smoothing parameter selection) in the
package `scam` from version 1.2-14. However, further careful development of
the stabilization strategies and convergence investigation are intended.

In practical terms, `scam` allows full fitting of exponential family GAMs with
shape restrictions on smooths, automatically estimating smoothing parameters
and uncertainty estimates. It provides a framework for generalized additive
modelling with a mixture of unconstrained terms and various shape-restricted
terms of both univariate and bivariate type (bivariate smooths can be subject
to shape constraints on either both or one of the covariates), which can also
accommodate parametric terms, simple random effects, varying coefficient
terms, the interaction of covariates, and also linear functionals with or
without shape constraints as model components.

## Acknowledgements

This work is partly supported by the Swedish Research Council (Reg. No.
2022-04190).

## References

  * Antoniadis et al. [2007] A. Antoniadis, J. Bigot, and I. Gijbels.  Penalized wavelet monotone regression.  _Statistics & probability letters_, 77(16):1608–1621, 2007.
  * Bollaerts et al. [2006] K. Bollaerts, P. Eilers, and I. van Mechelen.  Simple and multiple P-splines regression with shape constraints.  _British Journal of Mathematical and Statistical Psychology_ , 59:451–469, 2006.
  * Brezger and Steiner [2008] A. Brezger and W. J. Steiner.  Monotonic regression based on Bayesian P-splines: An application to estimating price response functions from store-level scanner data.  _Journal of business & economic statistics_, 26(1):90–104, 2008.
  * Chen and Samworth [2016] Y. Chen and R. J. Samworth.  Generalized additive and index models with shape constraints.  _JRSS (B)_ , 78(4):729–754, 2016.
  * Cheng [2009] G. Cheng.  Semiparametric additive isotonic regression.  _J Stat Plan Inference_ , 139(6):1980–1991, 2009\.
  * Cheng et al. [2012] G. Cheng, Y. Zhao, and B. Li.  Empirical likelihood inferences for the semiparametric additive isotonic regression.  _Journal of Multivariate Analysis_ , 112:172–182, 2012\.
  * Dette and Scheder [2006] H. Dette and R. Scheder.  Strictly monotone and smooth nonparametric regression for two or more variables.  _Canadian Journal of Statistics_ , 34(4):535–561, 2006.
  * Du et al. [2013] P. Du, C. F. Parmeter, and J. S. Racine.  Nonparametric kernel regression with multiple predictors and multiple shape constraints.  _Statistica Sinica_ , pages 1347–1371, 2013.
  * Dunson [2005] D. Dunson.  Bayesian semiparametric isotonic regression for count data.  _JASA_ , 100(470):618–627, 2005.
  * Dunson and Neelon [2003] D. Dunson and B. Neelon.  Bayesian inference on order-constrained parameters in glm.  _Biometrics_ , 59:286–295, 2003.
  * Eilers and Marx [1996] P. H. Eilers and B. D. Marx.  Flexible smoothing with b-splines and penalties.  _Statistical science_ , 11(2):89–121, 1996.
  * Fang and Meinshausen [2012] Z. Fang and N. Meinshausen.  Lasso isotone for high-dimensional additive isotonic regression.  _J Comput Graph Stat_ , 21(1):72–91, 2012.
  * Fellner [1986] W. H. Fellner.  Robust estimation of variance components.  _Technometrics_ , 28(1):51–60, 1986.
  * Groeneboom et al. [2001] P. Groeneboom, G. Jongbloed, and J. A. Wellner.  Estimation of a convex function: characterizations and asymptotic theory.  _The Annals of Statistics_ , 29(6):1653–1698, 2001.
  * Groeneboom et al. [2008] P. Groeneboom, G. Jongbloed, and J. A. Wellner.  The support reduction algorithm for computing non-parametric function estimates in mixture models.  _Scandinavian Journal of Statistics_ , 35(3):385–399, 2008.
  * Guntuboyina and Sen [2015] A. Guntuboyina and B. Sen.  Global risk bounds and adaptation in univariate convex regression.  _Probability Theory and Related Fields_ , 163(1):379–411, 2015.
  * Hall and Huang [2001] P. Hall and L.-S. Huang.  Nonparametric kernel regression subject to monotonicity constraints.  _The Annals of Statistics_ , 29(3):624–647, 2001\.
  * He and Shi [1998] X. He and P. Shi.  Monotone B-spline smoothing.  _JASA_ , 93(442):643–650, 1998.
  * Holmes and Heard [2003] C. Holmes and N. Heard.  Generalized monotonic regression using random change points.  _Statistics in Medicine_ , 22:623–638, 2003.
  * Kelly and Rice [1990] C. Kelly and J. Rice.  Monotone smoothing with application to dose-response curves and the assessment of synergism.  _Biometrics_ , 46:1071–1085, 1990.
  * Lang and Brezger [2004] S. Lang and A. Brezger.  Bayesian P-splines.  _J Comput Graph Stat_ , 13(1):183–212, 2004.
  * Lin and Dunson [2014] L. Lin and D. B. Dunson.  Bayesian monotone regression using gaussian process projection.  _Biometrika_ , 101(2):303–317, 2014.
  * Mammen [1991] E. Mammen.  Estimating a smooth monotone regression function.  _The Annals of Statistics_ , pages 724–740, 1991.
  * Mammen and Yu [2007] E. Mammen and K. Yu.  Additive isotone regression.  In _Asymptotics: particles, processes and inverse problems_ , pages 179–195. Institute of Mathematical Statistics, 2007.
  * Meyer [2018] M. C. Meyer.  A framework for estimation and inference in generalized additive models with shape and order restrictions.  _Statistical Science_ , 33(4):595–614, 2018.
  * Pya [2024] N. Pya.  _scam: Shape Constrained Additive Models_ , 2024.  URL https://CRAN.R-project.org/package=scam.  R package version 1.2-16.
  * Pya and Wood [2015] N. Pya and S. N. Wood.  Shape constrained additive models.  _Statistics and computing_ , 25:543–559, 2015.
  * Ramsay [1988] J. Ramsay.  Monotone regression splines in action.  _Stats. Science_ , 3(4):425–461, 1988.
  * Rousson [2008] V. Rousson.  Monotone fitting for developmental variables.  _Journal of Applied Statistics_ , 35(6):659–670, 2008.
  * Schall [1991] R. Schall.  Estimation in generalized linear models with random effects.  _Biometrika_ , 78(4):719–727, 1991.
  * Tutz and Leitenstorfer [2007] G. Tutz and F. Leitenstorfer.  Generalized smooth monotonic regression in additive modeling.  _J Comput Graph Stat_ , 16(1):165–188, 2007.
  * Villalobos and Wahba [1987] M. Villalobos and G. Wahba.  Inequality-constrained multivariate smoothing splines with application to the estimation of posterior probabilities.  _JASA_ , 82(397):239–248, 1987.
  * Wang and Meyer [2011] J. Wang and M. Meyer.  Testing the monotonicity or convexity of a function using regression splines.  _The Canadian Journal of Statistics_ , 39(1):89–107, 2011.
  * Wang and Xue [2015] L. Wang and L. Xue.  Constrained polynomial spline estimation of monotone additive models.  _J Stat Plan Inference_ , 167:27–40, 2015.
  * Wang and Zhou [2021] L. Wang and X.-H. Zhou.  Estimation of shape constrained additive models with missing response at random.  _Journal of Nonparametric Statistics_ , 33(1):118–133, 2021.
  * Wood [2023] S. Wood.  _mgcv: Mixed GAM Computation Vehicle with Automatic Smoothness Estimation_ , 2023.  URL https://CRAN.R-project.org/package=mgcv.  R package version 1.9.1.
  * Wood [2017] S. N. Wood.  _Generalized Additive Models: An Introduction with R_.  CRC press, 2017.
  * Wood [2020] S. N. Wood.  Inference and computation with generalized additive models and their extensions.  _Test_ , 29(2):307–339, 2020.
  * Wood and Fasiolo [2017] S. N. Wood and M. Fasiolo.  A generalized fellner-schall method for smoothing parameter optimization with application to tweedie location, scale and shape models.  _Biometrics_ , 73(4):1071–1081, 2017.
  * Wood et al. [2015] S. N. Wood, Y. Goude, and S. Shaw.  Generalized additive models for large data sets.  _JRSS (C)_ , 64(1):139–155, 2015.
  * Yu [2014] K. Yu.  On partial linear additive isotonic regression.  _Journal of the Korean Statistical Society_ , 43(1):11–17, 2014.
  * Zhang [2004] J. Zhang.  A simple and efficient monotone smoother using smoothing splines.  _Journal of Nonparametric Statistics_ , 16(5):779–796, 2004.