# High-dimensional Generalized Additive Mixed Model with Longitudinal Data

Mozhgan Taavoni  Department of Statistics, Faculty of Mathematical Sciences,
Ferdowsi University of Mashhad, IRAN<EMAIL_ADDRESS>Mohammad Arashi
Department of Statistics, Faculty of Mathematical Sciences, Ferdowsi
University of Mashhad, IRAN<EMAIL_ADDRESS>
###### Abstract

The problem of simultaneously estimating and selecting important variables is
challenging in the high-dimensional generalized additive mixed model (GAMM)
with longitudinal data because of the correlation structures. We give a
penalized estimation strategy for efficient variable selection when the number
of variables in the linear part and number of components in the non-linear
part can grow with the sample size. A Monte Carlo Newton-Raphson algorithm is
developed where a Metropolis step for generating random effects is used to
apply the proposed double-penalization technique in the high-dimensional GAMM.
In addition, to improve efficiency for regression coefficients, the estimation
of the working covariance matrix is involved in the proposed iterative
algorithm. We further study the asymptotic performance of the resulting
estimators and establish the oracle properties. We conduct extensive numerical
studies to assess the performance of our proposed estimation strategy and
numerically illustrate how efficient it is in selecting significant variables.

###### keywords:

Generalized estimating equations, Longitudinal data, Penalized likelihood,
Semiparametric additive mixed model, Variable selection, Oracle properties

##  1 Introduction

We often involve longitudinal studies in epidemiology, social science, and
other biomedical research areas, where repeated measurements from the same
subject are correlated. Under the likelihood principle, one may use the
generalized linear mixed models (GLMMs) to analyze longitudinal data (Zeger
and Karim, 1991; Breslow and Clayton, 1993), where subject-specific random
effects explicitly take care of the correlation induced by among-subject
variation. A standard GLMM assumes parametric fixed effects that may be too
restrictive when some functional covariate effects are present. Moreover, with
sparse data such as binary data, it is often not even possible to determine
empirically the underlying functional form. To allow richer and more flexible
model structures, an effective semiparametric regression tool is the
generalized additive mixed models (GAMMs) introduced by Lin and Zhang (1999),
which a linear predictor involves smooth functions of covariates. Statistical
inference of additive models with longitudinal data has also been considered
by some authors. This line of work includes Berhane and Tibshirani (1998), You
and Zhou (2007), Carroll et al. (2009), and Xue et al. (2010). However, when
the number of covariates is very large, selection of important variables is
still a challenging issue. Among all literature on variable selection we a can
refer to Li and Liang (2008) and Fan and Li (2012) where considered
semiparametric and mixed effects models, respectivilly. The literature on
variable selection for high-dimensional longitudinal data is rather limited
due to the challenges imposed by incorporating the intracluster correlation.
Some developments along this line can be found in Fan and Li (2004), Wang et
al. (2008), Ni et al. (2010), and Ma et al. (2013) for continuous longitudinal
data, and Fu (2003) and Dziak (2006) for discrete longitudinal data. In all
aforementioned studies, the dimension of predictors is fixed. Among few cases,
Wang et al. (2012), Chu et al. (2016), Chu et al. (2020), and Fang et al.
(2020) studied variable selection of the longitudinal GLMs with a diverging
number of parameters. The SCAD-penalized selection procedures were illustrated
in Xue et al. (2010) for the generalized additive model with correlated data.
To the best of our knowledge, simultaneous estimation and variable selection
in a high-dimensional GAMMs have not been investigated, especially when the
number of both linear predictors and functional components diverges with the
sample size. We try to fill this gap by allowing for non-Gaussian data and
nonlinear link functions.

###  1.1 Our contribution

Our model uses linear part and additive nonparametric functions to combine
underlying covariate effects while accounting for overdispersion and
correlation by adding random effects. We consider the case where the number of
linear predictors along with additive nonparametric functions are allowed to
increase with the sample size under a high-dimensional regime. We give a
regularization estimation procedure based on double penalization method while
the additive nonparametric functions are aproximated by smoothing splines, and
incorporate the within-cluster correlation to obtain an efficient parameter
estimation. We extend Taavoni and Arashi (2021), by applying penalty functions
to the estimating equation objective function to select the correct
generalized additive model. Our method is rather different from the latter
work because of including additive nonparametric components in the model and
using double penlization. Our numerical studies show that our proposal
performs well in estimation efficiency and model selection when the dimension
of linear predictors and functional components are high.

###  1.2 Plan of the paper

The rest of paper is organized as follows. Section 2, introduces the modeling
framework with the necessary notations, and spline approximation of the
nonparametric function. A double penalization procedure and the choice of the
tuning parameters included in Section 3. A classical SCAD penalty is used for
the linear part and the nonparametric functional forms are selected groupwise
for each functional component. In addition, the proposed penalized estimation
involves the specification of the posterior distribution of the random
effects, which cannot be evaluated in a closed form. However, it is possible
to approximate this posterior distribution by producing random draws from a
distribution using the Metropolis algorithm (Tanner, 1993), which does not
require the specification of the posterior distribution. In the same section,
we establish the asymptotic theory for the proposed method in a high-
dimensional framework, and achieve the oracle properties. To estimate the
parameters, a computationally flexible iterative algorithm is developed. In
Section 4, we apply a number of simulations to assess the finite sample
performance of the proposed estimation method in the GAMM. Real data analysis
is also presented in this section to augment the theoretical results. Some
concluding remarks are given in Section 5. Further, the proofs of the main
results as well as some instrumental lemmas are provided in a separate
supplementary file.

##  2 GAMMs and Approximation

Consider a longitudinal study with $n$ subjects and $n_{i}$ observations over
time for each subject with $i=1,\ldots,n$. Let $\boldsymbol{u}_{i}$ be a
$q\times 1$ vector of random effects corresponding to the $i$th subject,
$y_{ij}$, $j=1,\ldots,n_{i}$, denote the longitudinal response for subject $i$
measured at time $t_{ij}$,
$\boldsymbol{x}_{ij}=(x_{ij1},\ldots,x_{ijp_{n}})^{\top}$ is a $p_{n}\times 1$
vector of scalar covariates,
$\boldsymbol{z}_{ij}=(z_{ij1},\ldots,z_{ijq})^{\top}$ is a $q\times 1$ vector
of explanatory variables relating to the random effects. Suppose that
$\boldsymbol{y}_{i}=(y_{i1},\ldots,y_{in_{i}})^{\top}$ given
$\boldsymbol{u}_{i}$ are conditionally independent and each
$y_{ij}|\boldsymbol{b}_{i}$ is distributed according to an exponential
dispersion model with mean $\mu_{ij}=\mathbb{E}(y_{ij}|\boldsymbol{b}_{i})$
and conditional variance $\mathop{\rm
Var}\nolimits(y_{ij}|\boldsymbol{u}_{i})=\phi\omega^{-1}_{ij}\nu(\mu_{ij})$,
where $\nu(.)$ is a specified variance function, $\omega_{ij}$ is a known
weight (e.g., a binomial denominator), and $\phi$ is a dispersion parameter,
which can be either known or required to be estimated. In addition to random
and fixed effects, the model includes an additive term
$\boldsymbol{g}_{i}(x_{ijk})=\\{g_{i1}(x_{ij1}),\ldots,g_{ir_{n}}(x_{ijr_{n}})\\}$,
a $r_{n}\times 1$ vector of functional predictors. The $g_{ik}(x_{ijk})$ are
unknown smooth but arbitrary continuous and twice differentiable functions of
covariate $x_{ijk}$, on the bounded and closed interval $j_{g}$, with unknown
smooth mean functions $\mathbb{E}\big{\\{}g_{ik}(x_{ijk})\big{\\}}$ and
eigenvalue, eigenfunction pairs $\\{\pi_{kl},\phi_{kl}(x_{ijk})\\}_{l\geq 1}$,
$k=1,\ldots,r_{n}$. Then a GAMM is of the form

$\displaystyle
g(\mu_{ij})=\boldsymbol{x}_{ij}^{\top}\boldsymbol{\beta}_{n}+\sum_{k=1}^{r_{n}}g_{ik}(x_{ijk})+\boldsymbol{z}_{ij}^{\top}\boldsymbol{u}_{i},$
(1)

where $g(.)$ is a known monotonic and differential link function,
$\boldsymbol{x}_{ij}^{\top}\boldsymbol{\beta}_{n}$ is a linear parametric term
with parameter vector
$\boldsymbol{\beta}_{n}=(\beta_{1},\ldots,\beta_{p_{n}})^{\top}$,
$\sum_{k=1}^{r_{n}}g_{ik}(x_{ijk})$ is an additive term with unspecified
influence functions $g_{i1},\ldots,g_{ir_{n}}$ and finally
$\boldsymbol{z}_{ij}^{\top}\boldsymbol{b}_{i}$ contains the random effects
part. The $q$-dimensional random effects
$\boldsymbol{u}=\\{\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{q}\\}$ are
assumed to be independently and identically distributed, depending on
parameters $\boldsymbol{\Sigma}$ as $\boldsymbol{u}_{i}\sim
f_{u}(\boldsymbol{u}_{i}|\boldsymbol{\Sigma})$. If the functions $g_{k}(.)$
are strictly linear, the model reduces to the common GLMM. Note that the
dimension of random effects $q$ is assumed to be fixed, while the dimension of
both scalar covariates and functional predictors $p_{n}$ and $r_{n}$ may grow
to infinity with the sample size $n$. Specific assumptions will be considered
for the number of variables $p_{n}$ and $r_{n}$ in section 3.3.

In the estimation part, we approximate the smooth functions
$\\{g_{ik}(.)\\}_{k=1}^{r_{n}}$ by polynomial splines, which has been used by
Xue and Liang (2009), Xue et al. (2010), Groll and Tutz (2012), Wu and Xue
(2014) and Taavoni and Arashi (2021). For each $k=1,\ldots,p_{n}$, let
$\nu_{k}=\\{0=x_{k,0}<x_{k,1}<\ldots<x_{k,L_{nk}}<x_{k,L_{nk}+1}=1\\}$ be a
partition of the interval $[0,1]$, with $L_{nk}$ interior knots. The
polynomial splines of order $d+1$ are functions with $d$-degree (or less) of
polynomials on intervals $[x_{k,k^{{}^{\prime}}},x_{k,k^{{}^{\prime}}+1})$,
$k^{{}^{\prime}}=0,\ldots,L_{nk}-1$, and $[x_{k,L_{nk}},x_{k,L_{nk}+1}]$, and
have $p-1$ continuous derivatives globally. We denote the space of such spline
functions by $\phi_{k}=\big{(}\phi^{p}[0,1],\nu_{k}\big{)}$. Denote
$\phi_{k}^{0}=\\{g\in\phi_{k}:\int_{0}^{1}g(x)dx=0\\}$, which consists of
centered spline functions. Let $\\{B_{kl}(.)\\}_{l=1}^{h_{nk}}$ be a set of
spline bases of $\phi_{k}^{0}$ for variable $k$ with $h_{nk}=L_{nk}+d$. Here,
the centered truncated power basis is used in implementation, with
$\big{\\{}B_{kl}(.)=b_{kl}-\mathbb{E}_{n}(b_{kl})\big{\\}}_{l=1}^{h_{nk}}$,
where $\boldsymbol{b}_{k}=\\{b_{k1},\ldots,b_{kh_{nk}}\\}$ is the truncated
power basis given as
$\\{x_{k},\ldots,x_{k}^{d},(x_{k}-x_{k,1})_{+}^{d},\ldots,(x_{k}-x_{k,L_{nk}})_{+}^{d}\\}$,
in which $(x)_{+}=\max(0,x)$. Then each coefficient function $g_{ik}(.)$ for
$k=1,\ldots,r_{n}$ can be approximated by

$g_{ik}(x_{ijk})\thicksim\sum_{l=1}^{h_{nk}}\alpha_{kl}B_{kl}(x_{ijk}),$

where levels of is the spline coefficients. The subscript $k$ indicates that
different coefficient functions have different smoothness. Let
$\boldsymbol{\alpha}_{nk}=(\alpha_{k1},\ldots,\alpha_{kh_{nk}})^{\top}$ denote
the unknown parameter vector of the $j$th smooth function and
$\boldsymbol{B}_{k}(x_{ijk})=\big{(}B_{k1}(x_{ijk}),\ldots,B_{kh_{k}}(x_{ijk})\big{)}^{\top}$
represent the vector-valued evaluations of the $h_{nk}$ basis functions. Then,
model (1) has the parameterized form:

$\displaystyle
g(\mu_{ij})=\boldsymbol{x}_{ij}^{\top}\boldsymbol{\beta}_{n}+\boldsymbol{B}_{1}^{\top}(x_{ij1})\boldsymbol{\alpha}_{n1}+\ldots+\boldsymbol{B}_{r_{n}}^{\top}(x_{ijr_{n}})\boldsymbol{\alpha}_{nr_{n}}+\boldsymbol{z}_{ij}^{\top}\boldsymbol{u}_{i}.$

By collecting observations within one subject one obtains the design matrix
$\boldsymbol{X}_{i}=(\boldsymbol{x}_{i1},\ldots,\boldsymbol{x}_{in_{i}})^{\top}$
for the $i$th covariate, and analogously we set
$\boldsymbol{Z}_{i}=(\boldsymbol{z}_{i1},\ldots,\boldsymbol{z}_{in_{i}})^{\top}$,
so that the model has the simpler form

$\displaystyle
g(\boldsymbol{\mu}_{i})=\boldsymbol{X}_{i}\boldsymbol{\beta}_{n}+\boldsymbol{B}_{i1}\boldsymbol{\alpha}_{n1}+\ldots+\boldsymbol{B}_{ir_{n}}\boldsymbol{\alpha}_{nr_{n}}+\boldsymbol{Z}_{i}\boldsymbol{u}_{i},$

where
$\boldsymbol{B}_{ik}=\big{(}\boldsymbol{B}^{\top}_{1}(x_{i1k}),\ldots,\boldsymbol{B}^{\top}_{k}(x_{in_{i}k})\big{)}^{\top}$
denotes the B-spline design matrix of the $i$th subject and variable $k$ and
$g$ is understood componentwise. Furthermore, let
$\boldsymbol{X}=(\boldsymbol{X}_{1}^{\top},\ldots,\boldsymbol{X}_{n}^{\top})^{\top}$,
$\boldsymbol{Z}=\mathop{\rm
diag}\nolimits(\boldsymbol{Z}_{1},\ldots,\boldsymbol{Z}_{n})$ be a block-
diagonal matrix and,
$\boldsymbol{u}=(\boldsymbol{u}_{1}^{\top},\ldots,\boldsymbol{u}_{n}^{\top})^{\top}$
be the vector collecting all random effects. Then, it yields

$\displaystyle
g(\boldsymbol{\mu})=\boldsymbol{X}\boldsymbol{\beta}_{n}+\boldsymbol{B}_{1}\boldsymbol{\alpha}_{n1}+\ldots+\boldsymbol{B}_{r_{n}}\boldsymbol{\alpha}_{nr_{n}}+\boldsymbol{Z}\boldsymbol{u},$
(2)

where
$\boldsymbol{B}_{k}=(\boldsymbol{B}_{1k}^{\top},\ldots,\boldsymbol{B}_{nk}^{\top})^{\top}$
represens the B-spline design matrix of the $k$th smooth function. The model
can be more simplified to

$\displaystyle
g(\boldsymbol{\mu})=\boldsymbol{X}\boldsymbol{\beta}_{n}+\boldsymbol{B}\boldsymbol{\alpha}_{n}+\boldsymbol{Z}\boldsymbol{u},$

where
$\boldsymbol{\alpha}_{n}=(\boldsymbol{\alpha}_{n1}^{\top},\ldots,\boldsymbol{\alpha}_{nr_{n}}^{\top})^{\top}$
and
$\boldsymbol{B}=[\boldsymbol{B}_{1}^{\top},\ldots,\boldsymbol{B}_{r_{n}}^{\top}]^{\top}$.
An alternative form that we also use in the following is

$g(\boldsymbol{\mu})=\boldsymbol{D}\boldsymbol{\theta}_{n}+\boldsymbol{Z}\boldsymbol{u},$

where $\boldsymbol{D}=(\boldsymbol{X}^{\top},\boldsymbol{B}^{\top})^{\top}$
and
$\boldsymbol{\theta}_{n}=(\boldsymbol{\beta}_{n}^{\top},\boldsymbol{\alpha}_{n}^{\top})^{\top}$.

##  3 Regularization in the GAMMs

Focusing on the GLMM,we assume that the conditional density of
$y_{ij}|\boldsymbol{u}_{i}$, belongs to the exponential family whose
probability density function (p.d.f) is given by

$\displaystyle
p(y_{ij}|\boldsymbol{u}_{i})=\exp[\phi^{-1}\\{y_{ij}\theta_{ij}-b(\theta_{ij})\\}+c(y_{ij},\phi)],$

where $\theta_{ij}$ is the (scalar) canonical parameter, $b(\theta_{ij})$ is a
specific function corresponding to the type of exponential family, $c(.)$ the
log normalization constant and $\phi$ the dispersion parameter. Then the
quasi-likelihood (QL) of the model parameters
$(\boldsymbol{\theta}_{n},\boldsymbol{\Sigma},\phi)$ can be expressed as

$\displaystyle
L(\boldsymbol{\theta}_{n},\boldsymbol{\Sigma},\phi)=\prod_{i=1}^{n}\int
p_{\boldsymbol{y}_{i}|\boldsymbol{u}_{i}}(\boldsymbol{y}_{i}|\boldsymbol{u}_{i},\boldsymbol{\theta}_{n},\phi)p_{\boldsymbol{u}}(\boldsymbol{u}_{i}|\boldsymbol{\Sigma})d\boldsymbol{u}_{i},$
(3)

where
$p_{\boldsymbol{y}_{i}|\boldsymbol{u}_{i}}(\boldsymbol{y}_{i}|\boldsymbol{u}_{i},\boldsymbol{\theta}_{n},\phi)=\prod_{j=1}^{n_{i}}p(\boldsymbol{y}_{ij}|\boldsymbol{u}_{i},\boldsymbol{\theta}_{n},\phi)$.
In what follows, we consider the selection and estimation of both functional
and scalar parameters. Identifying important variables is a crucial step in
analyzing high-dimensional data, because each redundant variable involves an
infinite dimension of parameters for nonparametric components. Here a
predictor variable $x$ is said to be redundant in model (1), if and only if
$x=0$ or $g(x)=0$ almost surely. Otherwise, a predictor variable $x$ is said
to be relevant. Suppose that the true value $\boldsymbol{\beta}_{n0}$ of the
regression coefficient $\boldsymbol{\beta}_{n}$ can be decomposed into
$\boldsymbol{\beta}_{n0}=(\boldsymbol{\beta}_{n01}^{\top},\boldsymbol{\beta}_{n02}^{\top})^{\top}$
and the corresponding design matrix into
$\boldsymbol{X}_{i}=(\boldsymbol{X}_{i(1)},\boldsymbol{X}_{i(2)})$, where
$\boldsymbol{\beta}_{n01}$ is a $q_{n}\times 1$ vector corresponding to
significant covariate effects $\boldsymbol{X}_{i(1)}$ and
$\boldsymbol{\beta}_{n02}$ is a $(p_{n}-q_{n})\times 1$ vector of zeros. Also
assume that only the first $d_{n}$ functional predictors are significant or,
equivalently that the true values of $g_{ik}(x_{ijk})$ are such that
$g_{ik}(x_{ijk})\equiv 0$ for $k=d_{n}+1,\ldots,r_{n};x_{ijk}\in J_{x}$. In
other words,
$\boldsymbol{\alpha}_{n(d_{n}+1)},\ldots,\boldsymbol{\alpha}_{nr_{n}}$ are
zero vectors. Our goal is to consistently identify such subsets of relevant
variables and estimate their unknown parameters or function components. In
order to address the difficulty caused by the infinite dimensionality of
functional predictors and high dimensionality of scalar covariates, the terms
$g_{ik}(x_{ijk})$ in 1 are first approximated with truncated linear predictors
$\sum_{l=1}^{h_{nk}}\alpha_{kl}B_{kl}(x_{ijk})$ for some large numbers
$\\{h_{nk}\\}_{k=1}^{r_{n}}$. The numbers $h_{nk}$ may vary with the sample
size $n$, and play the role of smoothing parameters which balance the trade-
off between bias and variance. Then, a suitable regularization is imposed on
functional predictor $\boldsymbol{g}_{i}(.)$ and scalar covariate
$\boldsymbol{X}$, respectively. As a consequence, the coefficient functions of
unimportant functional predictors and regression parameters of unimportant
scalar covariates will be shrunken to zero. Therefore, the goals of the
variable selection and estimation of important functional coefficients and
fixed-effects regressor coefficients can be achieved by maximizing the
penalized QL

$\displaystyle\ell(\boldsymbol{\theta}_{n},\boldsymbol{\Sigma},\phi)$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\hbox{ln}p_{\boldsymbol{y}_{i}|\boldsymbol{u}_{i}}(\boldsymbol{y}_{i}|\boldsymbol{u}_{i},\boldsymbol{\theta}_{n},\phi)+\sum_{i=1}^{n}\hbox{ln}p_{\boldsymbol{u}_{i}}(\boldsymbol{u}_{i}|\boldsymbol{\Sigma})$
$\displaystyle-$ $\displaystyle
n\sum_{k=1}^{p_{n}}p_{\lambda_{n}}(|\beta_{nk}|)-n\sum_{k=1}^{p_{n}}p_{\lambda_{nk}}(\|\boldsymbol{\alpha}_{nk}\|_{\boldsymbol{W}_{k}}),$

where $p_{\lambda_{n}}(.)$ and $p_{\lambda_{nk}}(.)$ are penalty functions
with tuning parameters $\lambda_{n}$, $\lambda_{nk}$, and
$\|\boldsymbol{\alpha}_{nk}\|_{\boldsymbol{W}_{k}}=\sqrt{\boldsymbol{\alpha}_{nk}^{\top}\boldsymbol{W}_{k}\boldsymbol{\alpha}_{nk}}$
with
$\boldsymbol{W}_{k}=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{n_{i}}\sum_{j=1}^{n_{i}}\boldsymbol{B}_{k}(\boldsymbol{x}_{ijk})\boldsymbol{B}^{\top}_{k}(\boldsymbol{x}_{ijk})$,
the weighted Euclidean norm associated with a group penalty. In order to
regularize all predictors on a comparable scale, one often standardizes the
predictors before imposing a penalty on the objective function. Thus, we
standardize $(\boldsymbol{x}_{1k},\ldots,\boldsymbol{x}_{nk})$ to have unit
variance, for $k=1,\ldots,p_{n}$. Since the variability of the $k$th
functional predictor $g_{ik}(x_{ijk})$ can be approximated by
$\sum_{l=1}^{h_{nk}}\widehat{\pi}_{kl}$, and standardization of
$g_{ik}(x_{ijk})$ is equivalent to adding weight to the penalty function, we
suggest using
$p_{\lambda_{kn}}=\lambda_{n}\Big{(}\sum_{l=1}^{h_{nk}}\widehat{\pi}_{kl}\Big{)}^{1/2}$,
which simplifies both the computation and theory. The resulting estimators for
$\beta_{nk}$ and $\boldsymbol{\alpha}_{nk}$ obtained from (3) are then denoted
by $\widehat{\beta}_{nk}$ and $\widehat{\boldsymbol{\alpha}}_{nk}$. For the
maximum QL (MQL) estimates of the parameters
$(\boldsymbol{\theta}_{n},\boldsymbol{\Sigma},\phi)$, one can maximize the
function (3) by using suitable numerical techniques.

The EM algorithm is an attractive method to obtain the ML estimates, in
presence of incomplete data, which avoids explicit calculation of the observed
data log-likelihood. To set up the EM algorithm we consider the random
effects, $\boldsymbol{u}_{i}$, to be the missing data and the complete data,
is then $(\boldsymbol{y}_{i},\boldsymbol{u}_{i})$. Considering
$\boldsymbol{u}_{i}$ to be the missing has an advantage that in the M-step,
maximization can be accomplished with respect to the parameters
$\boldsymbol{\theta}_{n}$ and $\phi$ in the first, third, and fourth terms of
(3). Thus, the M-step with respect to $(\boldsymbol{\theta}_{n},\phi)$ uses
only the GLM part of the penalized likelihood function. Hence, the procedure
is similar to a standard GLM computation assuming $\boldsymbol{u}_{i}$ is
known. Therefore, maximizing with respect to $\boldsymbol{\Sigma}$, in the
second term, can be handled by the ML using the distribution of
$p_{\boldsymbol{u}_{i}}(\boldsymbol{u}_{i}|\boldsymbol{\Sigma})$ after
replacing sufficient statistics with the conditional expected values. It thus
makes sense to develop a simulation analogous to the Newton-Raphson approach
for fitting the GAMM. Using this separation as in (3) and the fact that
$\frac{\partial
p_{\lambda}(|\beta_{nk}|)}{\partial\beta_{nk}}=p^{{}^{\prime}}_{\lambda}(|\beta_{nk}|)\hbox{sign}(\beta_{nk})$
and $\frac{\partial
p_{\lambda}(\|\boldsymbol{\alpha}_{nk}\|_{\boldsymbol{W}_{k}})}{\partial\boldsymbol{\alpha}_{nk}}=p^{{}^{\prime}}_{\lambda}(\|\boldsymbol{\alpha}_{nk}\|_{\boldsymbol{W}_{k}})\|\boldsymbol{\alpha}_{nk}\|_{\boldsymbol{W}_{k}}^{-1}\boldsymbol{W}_{k}\boldsymbol{\alpha}_{nk}$,
the penalized MQL (PMQL) equations for $\boldsymbol{\theta}_{n}$ and
$\boldsymbol{\Sigma}$ take the following forms

$\displaystyle\mathbb{E}\Big{[}\frac{\partial\ln
p_{\boldsymbol{y}_{i}|\boldsymbol{u}_{i}}(\boldsymbol{y}_{i}|\boldsymbol{\theta}_{n},\phi)}{\partial\boldsymbol{\theta}_{n}}\Big{|}y_{ij}\Big{]}-\sum_{k=1}^{p_{n}}q_{\lambda_{n}}(|\boldsymbol{\beta}_{n}|)\hbox{sign}(\boldsymbol{\beta}_{n})$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad-\sum_{k=1}^{p_{n}}q_{\lambda_{nk}}(\|\boldsymbol{\alpha}_{nk}\|_{\boldsymbol{W}_{k}})\|\boldsymbol{\alpha}_{nk}\|_{\boldsymbol{W}_{k}}^{-1}\boldsymbol{W}_{k}\boldsymbol{\alpha}_{nk}=\boldsymbol{0};$
$\displaystyle\mathbb{E}\Big{[}\frac{\partial\ln
p_{\boldsymbol{u}_{i}}(\boldsymbol{u}_{i}|\boldsymbol{\Sigma})}{\partial\boldsymbol{\Sigma}}|y_{ij}\Big{]}=\boldsymbol{0},$

with $\hbox{sign}(a)=I(a>0)-I(a<0)$ and $q_{\lambda}(.)$ is the first-order
derivative of $p_{\lambda}(.)$. The PML estimates of
$(\boldsymbol{\beta}_{n},\boldsymbol{\alpha}_{n})$ and $\boldsymbol{\Sigma}$
can be obtained by solving the preceding equations numerically. Following
Taavoni and Arashi (2021), the optimal penalized generalized estimating
equation (PGEE) for $\boldsymbol{\theta}_{n}$ is given by

$\displaystyle\mathbb{E}_{\boldsymbol{u}|\boldsymbol{y}}\Big{[}n^{-1}\sum_{i=1}^{n}\frac{\partial\boldsymbol{\mu}_{i}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})}{\partial\boldsymbol{\theta}_{n}^{\top}}\boldsymbol{V}_{i}^{-1}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\big{(}\boldsymbol{y}_{i}-\boldsymbol{\mu}_{i}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\big{)}\Big{]}$
(5)
$\displaystyle-\sum_{k=1}^{p_{n}}q_{\lambda_{n}}(|\boldsymbol{\beta}_{n}|)\hbox{sign}(\boldsymbol{\beta}_{n})-\sum_{k=1}^{p_{n}}q_{\lambda_{nk}}(\|\boldsymbol{\alpha}_{nk}\|_{\boldsymbol{W}_{k}})\|\boldsymbol{\alpha}_{nk}\|_{\boldsymbol{W}_{k}}^{-1}\boldsymbol{W}_{k}\boldsymbol{\alpha}_{nk}=\boldsymbol{0},$

where
$\boldsymbol{\mu}_{i}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})=(\mu_{i1},\ldots,\mu_{in_{i}})^{\top}$
and $\boldsymbol{V}_{i}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})$ is the
covariance matrix of $\boldsymbol{y}_{i}|\boldsymbol{u}_{i}$. In real
applications the true intracluster covariance structure is often unknown. We
assume the working correlation matrix
$\boldsymbol{R}(\boldsymbol{\rho}):\boldsymbol{V}_{i}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})=\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\boldsymbol{R}(\boldsymbol{\rho})\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})$,
where $\boldsymbol{\rho}$ is a finite dimensional parameter and
$\boldsymbol{A}_{i}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})=\mathop{\rm
diag}\nolimits(\nu_{i1},\ldots,\nu_{in_{i}})$. Some commonly used working
correlation structures include independence, autocorrelation (AR)-1,
exchangeable, toeplitz, M-dependent, or unstructured correlation, among
others. For a given working correlation structure, $\boldsymbol{\rho}$ can be
estimated using the residual-based method of moments. With the estimated
working correlation matrix
$\widehat{\boldsymbol{R}}\equiv\boldsymbol{R}(\widehat{\boldsymbol{\rho}})$,
the PGEE in (5) reduces to

$\displaystyle\boldsymbol{U}_{n}(\boldsymbol{\theta}_{n})=\boldsymbol{S}_{n}(\boldsymbol{\theta}_{n})-\sum_{k=1}^{p_{n}}q_{\lambda_{n}}(|\boldsymbol{\beta}_{n}|)\hbox{sign}(\boldsymbol{\beta}_{n})-\sum_{k=1}^{p_{n}}q_{\lambda_{nk}}(\|\boldsymbol{\alpha}_{nk}\|_{\boldsymbol{W}_{k}})\|\boldsymbol{\alpha}_{nk}\|_{\boldsymbol{W}_{k}}^{-1}\boldsymbol{W}_{k}\boldsymbol{\alpha}_{nk}=\boldsymbol{0},$

where
$\boldsymbol{S}_{n}(\boldsymbol{\theta}_{n})=\mathbb{E}_{u|y}\Big{[}\sum_{i=1}^{n}\boldsymbol{D}_{i}^{\top}\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\widehat{\boldsymbol{R}}^{-1}\boldsymbol{A}_{i}^{-\frac{1}{2}}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\big{(}\boldsymbol{y}_{i}-\boldsymbol{\mu}_{i}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\big{)}\Big{]}$.
We formally define the estimator as the solution
$\widehat{\boldsymbol{\theta}}_{n}=(\widehat{\boldsymbol{\beta}}^{\top},\widehat{\boldsymbol{\alpha}}^{\top})^{\top}$
of the above estimating equations.

Among all penalty functions, the smoothing clipped absolute deviation (SCAD)
penalty proposed by Fan and Li (2001) can be used to retain the good features
of both subset selection and ridge regression, for producing sparse solutions,
and to ensure continuity of the selected models. Therefore, we will use the
SCAD penalty in our simulation and application studies. The SCAD penalty
function is defined by

$q_{\lambda}(\theta)=p^{{}^{\prime}}_{\lambda_{n}}(\theta)=\lambda\Big{\\{}I(\theta\leq\lambda)+\frac{(a\lambda-\theta)_{+}}{(a-1)\lambda}I(\theta>\lambda)\Big{\\}};\quad
a>2,$

where the notation $(.)_{+}$ stands for the positive part of $(.)$.

###  3.1 Computational algorithm

Using the local quadratic approximation (LQA), in the neighborhoods of the
true parameter values $\beta_{n0k}$, $|\beta_{n0k}|>0$, the derivative of the
penalty function is well approximated by
$q_{\lambda_{n}}(|\beta_{nk}|)\hbox{sign}(\beta_{nk})\thickapprox\frac{q_{\lambda_{n}}(|\beta_{n0k}|)}{|\beta_{n0k}|}\beta_{nk}$.
If $\boldsymbol{\alpha}_{nk}$ is close to $0$ in the sense that
$\|\boldsymbol{\alpha}_{nk}\|_{\boldsymbol{W}_{k}}\leq\epsilon$, for some
small threshold value $\epsilon$, set $\boldsymbol{\alpha}_{nk}$ to
$\boldsymbol{0}$. In the implementation, we have used $\epsilon=10^{-6}$. Xue
(2009) and Xue et al. (2010), applied the following LQA for the nonzero
functional component $\boldsymbol{\alpha}_{nk}$, that for an initial value
$\|\boldsymbol{\alpha}_{n0k}\|_{\boldsymbol{W}_{k}}>\epsilon$, one can locally
approximate the penalty function by

$\displaystyle
p_{\lambda_{nk}}(\|\boldsymbol{\alpha}_{nk}\|_{\boldsymbol{W}_{k}})$
$\displaystyle\thickapprox$ $\displaystyle
p_{\lambda_{nk}}(\|\boldsymbol{\alpha}_{n0k}\|_{\boldsymbol{W}_{k}})$
$\displaystyle+$
$\displaystyle\frac{1}{2}p^{{}^{\prime}}_{\lambda_{nk}}(\|\boldsymbol{\alpha}_{n0k}\|_{\boldsymbol{W}_{k}})\|\boldsymbol{\alpha}_{n0k}\|^{-1}_{\boldsymbol{W}_{k}}(\boldsymbol{\alpha}^{\top}_{nk}\boldsymbol{W}_{k}\boldsymbol{\alpha}_{nk}-\boldsymbol{\alpha}^{\top}_{n0k}\boldsymbol{W}_{k}\boldsymbol{\alpha}_{n0k}).$

With the mentioned LQAs, we apply the Newton–Raphson method to solve
$\boldsymbol{U}_{n}(\widehat{\boldsymbol{\theta}}_{n})=0$, and get the
following updating formula

$\displaystyle\widehat{\boldsymbol{\theta}}_{n}^{(m+1)}=\widehat{\boldsymbol{\theta}}_{n}^{(m)}+\Big{\\{}\boldsymbol{H}_{n}(\widehat{\boldsymbol{\theta}}_{n}^{(m)})+n\boldsymbol{E}_{n}(\widehat{\boldsymbol{\theta}}_{n}^{(m)})\Big{\\}}^{-1}\times\Big{\\{}\boldsymbol{S}_{n}(\widehat{\boldsymbol{\theta}}_{n}^{(m)})+n\boldsymbol{E}_{n}(\widehat{\boldsymbol{\theta}}_{n}^{(m)})\widehat{\boldsymbol{\theta}}_{n}^{(m)}\Big{\\}},$
(6)

where

$\displaystyle\boldsymbol{H}_{n}(\widehat{\boldsymbol{\theta}}_{n}^{(m)})=\mathbb{E}_{u|y}\Big{[}\sum_{i=1}^{n}\boldsymbol{D}_{i}^{\top}\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\widehat{\boldsymbol{R}}^{-1}\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\boldsymbol{D}_{i}\Big{]},$
$\displaystyle\boldsymbol{E}_{n}(\widehat{\boldsymbol{\theta}}_{n}^{(m)})=\mathop{\rm
diag}\nolimits\Big{\\{}\frac{q_{\lambda_{n}}(|\beta_{n1}|)}{\epsilon+|\beta_{n1}|},\ldots,\frac{q_{\lambda_{n}}(|\beta_{np_{n}}|)}{\epsilon+|\beta_{np_{n}}|},$
$\displaystyle
q_{\lambda_{n1}}(\|\boldsymbol{\alpha}_{n1}\|_{\boldsymbol{W}_{1}})\|\boldsymbol{\alpha}_{n1}\|_{\boldsymbol{W}_{1}}^{-1},\ldots,$
$\displaystyle
q_{\lambda_{nr_{n}}}(\|\boldsymbol{\alpha}_{nr_{n}}\|_{\boldsymbol{W}_{r_{n}}})\|\boldsymbol{\alpha}_{nr_{n}}\|_{\boldsymbol{W}_{r_{n}}}^{-1}\Big{\\}}.$

After obtaining the estimator
$\widehat{\boldsymbol{\theta}}_{n}=(\widehat{\boldsymbol{\beta}}_{n}^{\top},\widehat{\boldsymbol{\alpha}}_{n}^{\top})^{\top}$
through penalization in (6), for any given $\boldsymbol{x}\in[0,1]^{r_{n}}$,
an estimator of the unknow functional components in (1) is given as

$\displaystyle\widehat{g}_{k}(x_{ijk})=\sum_{l=1}^{h_{nk}}\widehat{\boldsymbol{\alpha}}_{kl}B_{kl}(x_{ijk}).$

Note that, in general, the expectations in (6) cannot be computed in a closed
form as the conditional distribution of
$\boldsymbol{u}_{i}|\boldsymbol{y}_{i}$ involves the marginal distribution of
$\boldsymbol{y}_{i}$, which is not easy to be computed explicitly. In the
following we outline the computational procedure used for sample generation.
Let $\boldsymbol{U}$ denote the previous draw from the conditional
distribution of $\boldsymbol{U}|\boldsymbol{y}$, and generate a new value
$u_{k}^{*}$ for the $j$th component of
$\boldsymbol{U}^{*}=(u_{1},\ldots,u_{k-1},u^{*}_{k},u_{k+1},\ldots,u_{nq})$ by
using the candidate distribution $p_{\boldsymbol{u}}$, accept
$\boldsymbol{U}^{*}$ as the new value with probability

$\displaystyle\alpha_{k}(\boldsymbol{U},\boldsymbol{U}_{*})=\hbox{min}\Big{\\{}1,\frac{p_{\boldsymbol{u}|\boldsymbol{y}}(\boldsymbol{U}^{*}|\boldsymbol{y},\boldsymbol{\theta}_{n},\boldsymbol{D})p_{u}(\boldsymbol{U}|\boldsymbol{D})}{p_{\boldsymbol{u}|\boldsymbol{y}}(\boldsymbol{U}|\boldsymbol{y},\boldsymbol{\theta}_{n},\boldsymbol{D})p_{u}(\boldsymbol{U}^{*}|\boldsymbol{D})}\Big{\\}}.$
(7)

otherwise, reject it and retain the previous value $\boldsymbol{U}$. The
second term in brace in (7) can be simplified to

$\displaystyle\frac{p_{\boldsymbol{u}|\boldsymbol{y}}(\boldsymbol{U}^{*}|\boldsymbol{y},\boldsymbol{\theta}_{n},\boldsymbol{D})p_{u}(\boldsymbol{U}|\boldsymbol{D})}{p_{\boldsymbol{u}|\boldsymbol{y}}(\boldsymbol{U}|\boldsymbol{y},\boldsymbol{\theta}_{n},\boldsymbol{D})p_{u}(\boldsymbol{U}^{*}|\boldsymbol{D})}$
$\displaystyle=$
$\displaystyle\frac{p_{\boldsymbol{y}|\boldsymbol{u}}(\boldsymbol{y}|\boldsymbol{U}^{*},\boldsymbol{\theta}_{n})}{f_{\boldsymbol{y}|\boldsymbol{u}}(\boldsymbol{y}|\boldsymbol{U},\boldsymbol{\theta}_{n})}$
$\displaystyle=$
$\displaystyle\frac{\prod_{i=1}^{n}p_{\boldsymbol{y}_{i}|\boldsymbol{u}}(\boldsymbol{y}_{i}|\boldsymbol{U}^{*},\boldsymbol{\theta}_{n})}{\prod_{i=1}^{n}f_{\boldsymbol{y}_{i}|\boldsymbol{u}}(\boldsymbol{y}_{i}|\boldsymbol{U},\boldsymbol{\theta}_{n})}.$

Note that, the calculation of the acceptance function
$\alpha_{k}(\boldsymbol{U},\boldsymbol{U}_{*})$ here involves only the
specification of the conditional distribution of
$\boldsymbol{y}|\boldsymbol{u}$ which can be computed in a closed form.

The suggested double penalized approach can be implemented step by step, and
the detailed algorithm can be summarized in Algorithm 1 descring the
Metropolis step into the Newton-Raphson iterative equation (6) for the Monte
Carlo estimates of expected values.

Algorithm 1 Monte Carlo Newton-Raphson (MCNR) algorithm

step 1. Smooth each predictor trajectory $\\{g_{ik}(.)\\}_{k=1}^{p_{n}}$, by
local linear smoothing technique. The smoothed predictor trajectories are then
denoted as $g_{ik}(.)=\boldsymbol{B}_{ik}\boldsymbol{\alpha}_{nk}$. Set
$m_{k}=0$. Choose initial values
$\boldsymbol{\theta}_{n}^{0}=(\boldsymbol{\beta}_{n}^{0\top},\boldsymbol{\alpha}_{n}^{0\top})^{\top}$
and $\boldsymbol{\Sigma}^{0}$.

step 2. Generate $N$ observations
$\boldsymbol{U}^{(1)},\ldots,\boldsymbol{U}^{(N)}$ from the distribution
$p_{\boldsymbol{u}|\boldsymbol{y}}(\boldsymbol{u}|\boldsymbol{y},\boldsymbol{\theta}_{n}^{(m_{k})},\boldsymbol{\Sigma}^{(m_{k})})$
using the Metropolis algorithm. Use these observations to find the Monte Carlo
estimates of the expectations. Specially,

  * a)

Compute $\boldsymbol{\theta}_{n}^{(m_{k}+1)}$ from the expression

$\displaystyle\boldsymbol{\theta}_{n}^{(m_{k}+1)}$ $\displaystyle=$
$\displaystyle\boldsymbol{\theta}_{n}^{(m_{k})}+\Big{\\{}\frac{1}{N}\sum_{k=1}^{N}\Big{[}\boldsymbol{H}_{n}\big{(}\widehat{\boldsymbol{\theta}}_{n}^{(m_{k})},\boldsymbol{U}^{(k)}\big{)}\Big{]}+n\boldsymbol{E}_{n}(\widehat{\boldsymbol{\theta}}_{n}^{(m_{k})})\Big{\\}}^{-1}$
$\displaystyle\qquad\times\Big{\\{}\frac{1}{N}\sum_{k=1}^{N}\Big{[}\boldsymbol{S}_{n}(\widehat{\boldsymbol{\theta}}_{n}^{(m_{k})},\boldsymbol{U}^{(k)})\Big{]}-n\boldsymbol{E}_{n}(\widehat{\boldsymbol{\theta}}_{n}^{(m_{k})})\widehat{\boldsymbol{\theta}}_{n}^{(m_{k})}\Big{\\}},$

where

$\displaystyle\boldsymbol{H}_{n}\big{(}\widehat{\boldsymbol{\theta}}_{n}^{(m_{k})},\boldsymbol{U}^{(k)}\big{)}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\boldsymbol{D}_{i}^{\top}\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n}^{(m_{k})},U_{i}^{(k)})\widehat{\boldsymbol{R}}^{-1}\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n}^{(m_{k})},U_{i}^{(k)})\boldsymbol{D}_{i},$
$\displaystyle\boldsymbol{S}_{n}\big{(}\widehat{\boldsymbol{\theta}}_{n}^{(m_{k})},\boldsymbol{U}^{(k)}\big{)}$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{n}\boldsymbol{D}_{i}^{\top}\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n}^{(m_{k})},U_{i}^{(k)})\widehat{\boldsymbol{R}}^{-1}\boldsymbol{A}_{i}^{-\frac{1}{2}}(\boldsymbol{\theta}_{n}^{(m_{k})},U_{i}^{(k)})$
$\displaystyle\qquad\qquad\qquad\times\big{(}\boldsymbol{y}_{i}-\boldsymbol{\mu}_{i}(\boldsymbol{\theta}_{n}^{(m_{k})},U_{i}^{(k)})\big{)}.$

  * b)

Compute $\boldsymbol{\Sigma}^{(m_{k}+1)}$ by maximizing
$\frac{1}{N}\sum_{k=1}^{N}\hbox{ln}f_{u}(\boldsymbol{U}^{(k)}|\boldsymbol{\Sigma}).$

  * c)

Set $m_{k}=m_{k}+1$.

step 3. Go to step 2 until convergence is achieved. Choose
$\boldsymbol{\theta}_{n}^{(m_{k}+1)}$ and $\boldsymbol{\Sigma}^{(m_{k}+1)}$ to
be the MCNR estimates of $\boldsymbol{\theta}_{n}$ and $\boldsymbol{\Sigma}$.

To implement the proposed double penalized procedure (5), one needs to choose
appropriate spline spaces $\\{\phi_{k}^{n}\\}_{k=1}^{p_{n}}$ and tuning
parameters $a$, $\lambda_{n}$ involved in the SCAD penalty. For the choices of
$\\{\phi_{k}^{n}\\}_{k=1}^{p_{n}}$ we use splines with equally spaced knots
and fixed degrees, and select only $L_{nk}$, the number of interior knots
using the data. For computational convenience, we use equally spaced knots
with the number of interior knots $L_{n}\approx n^{1/(2r+1)}$ following Xue
(2009) and Xue et al. (2010). Set $a=3.7$ and finally, choose $\lambda_{n}$.
The tuning parameter $\lambda_{n}$ controls the sparsity of both the
functional and the regression coefficients, then to obtain the optimal value
of $\lambda_{n}$, one can modify the generalized cross validation (GCV)
procedure. Due to the lack of joint likelihood in the generalized model, to
select the tuning parameter $\lambda_{n}$, we use GCV suggested by Fan and Li
(2001) as defined by

$\displaystyle\hbox{GCV}_{\lambda_{n}}=\frac{\hbox{RSS}(\lambda_{n})/n}{(1-d(\lambda_{n})/n)^{2}},$

where

$\displaystyle\hbox{RSS}(\lambda_{n})=\frac{1}{N}\sum_{k=1}^{N}\Big{[}\sum_{i=1}^{n}\big{(}\boldsymbol{y}_{i}-\boldsymbol{\mu}_{i}(\widehat{\boldsymbol{\theta}}_{n},U_{i}^{(k)})\big{)}^{\top}\boldsymbol{W}_{i}^{-1}\big{(}\boldsymbol{y}_{i}-\boldsymbol{\mu}_{i}(\widehat{\boldsymbol{\theta}}_{n},U_{i}^{(k)})\big{)}\Big{]}$
(8)

is the residual sum of squares, and

$\displaystyle
d(\lambda_{n})=tr\Big{[}\Big{\\{}\frac{1}{N}\sum_{k=1}^{N}\Big{[}\boldsymbol{H}_{n}\big{(}\widehat{\boldsymbol{\theta}}_{n},\boldsymbol{U}^{(k)}\big{)}\Big{]}+n\boldsymbol{E}_{n}(\widehat{\boldsymbol{\theta}}_{n})\Big{\\}}^{-1}\times\Big{\\{}\frac{1}{N}\sum_{k=1}^{N}\Big{[}\boldsymbol{H}_{n}\big{(}\widehat{\boldsymbol{\theta}}_{n},\boldsymbol{U}^{(k)}\big{)}\Big{]})\Big{\\}}\Big{]}$

is the effective number of parameters. Then, $\lambda_{opt}$ is the minimizer
of the $\hbox{GCV}_{\lambda_{n}}$. Note that $\boldsymbol{W}_{i}$ in (8) is an
$n_{i}\times n_{i}$ covariance matrix of $\boldsymbol{y}_{i}$, that can be
computed as
$\boldsymbol{W}_{i}=\mathbb{E}_{u|y}\Big{(}\hbox{var}(\boldsymbol{y}_{i}|\boldsymbol{u}_{i})\Big{)}+\hbox{var}_{u|y}\Big{(}\mathbb{E}(\boldsymbol{y}_{i}|\boldsymbol{u}_{i})\Big{)},$
where

$\displaystyle\mathbb{E}_{u|y}\Big{(}\hbox{var}(\boldsymbol{y}_{i}|\boldsymbol{u}_{i})\Big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{k=1}^{N}\Big{[}\boldsymbol{V}_{i}(\widehat{\boldsymbol{\theta}}_{n},U_{i}^{(k)})\Big{]},$
$\displaystyle\hbox{var}_{u|y}\Big{(}\mathbb{E}(\boldsymbol{y}_{i}|\boldsymbol{u}_{i})\Big{)}$
$\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{k=1}^{N}\Big{[}\boldsymbol{\mu}_{i}(\widehat{\boldsymbol{\theta}}_{n},U_{i}^{(k)})\Big{]}^{2}-\bigg{[}\frac{1}{N}\sum_{k=1}^{N}\Big{[}\boldsymbol{\mu}_{i}(\widehat{\boldsymbol{\theta}}_{n},U_{i}^{(k)})\Big{]}\bigg{]}^{2}.$

###  3.2 Oracle properties

Our proposed estimator for $\boldsymbol{\theta}_{n}$ is the solution of
$\boldsymbol{U}_{n}(\boldsymbol{\theta}_{n})=\boldsymbol{0}$. Because
$\boldsymbol{U}_{n}(\boldsymbol{\theta}_{n})$ has discontinuous points, an
exact solution to $\boldsymbol{U}_{n}(\boldsymbol{\theta}_{n})=\boldsymbol{0}$
may not exist. We formally define the estimator
$\widehat{\boldsymbol{\theta}}_{n}$ to be an approximate solution, i.e.,
$\boldsymbol{U}_{n}(\widehat{\boldsymbol{\theta}}_{n})=o(a_{n})$ for a
sequence $a_{n}\to 0$. The rate of $a_{n}$ will be specified in Theorem 1.
Meanwhile, If Eq. (5) has multiple solutions, only a sequence of consistent
estimator $\widehat{\boldsymbol{\theta}}_{n}$ is considered. A sequence
$\widehat{\boldsymbol{\theta}}_{n}$ is said to be a consistent sequence, if
$\widehat{\boldsymbol{\beta}}_{n}-\boldsymbol{\beta}_{n0}\to\boldsymbol{0}$
and
$\underset{x}{\sup}|\boldsymbol{B}^{\top}(x)\widehat{\boldsymbol{\alpha}}_{n}-g_{0}(x)|\to\boldsymbol{0}$
in probability as $n\to\infty$. In what follows, we establish the asymptotic
properties of the estimators of the parametric and nonparametric components
for the GAMM under the proposed PGEE. In our asymptotic study, we assume the
number of subjects $n$ goes to infinity, $\phi=1$ and $n_{i}=m<\infty$.
Extension of the methodology to the cases of unequal $n_{i}$ is
straightforward. We vary the dimension of $\boldsymbol{A}_{i}$ and replace
$\widehat{\boldsymbol{R}}$ by $\widehat{\boldsymbol{R}}_{i}$, which is the
$n_{i}\times n_{i}$ matrix using the specified working correlation structure
and the corresponding initial parameter $\boldsymbol{\rho}$ estimator. For
simplicity, we further assume that the response, the functional, random and
scalar predictors have been centered to mean zero for a total of
$N=\sum_{i=1}^{n_{i}}$ observation. Without loss of generality, we denote
$g(\mu_{ij})=\boldsymbol{X}_{i(1)}\boldsymbol{\beta}_{n1}+\boldsymbol{X}_{i(2)}\boldsymbol{\beta}_{n2}+\sum_{k=1}^{s_{n}}g_{ik}(x_{ijk})+\sum_{k=s_{n}+1}^{p_{n}}g_{ik}(x_{ijk})+\boldsymbol{z}_{ij}^{\top}\boldsymbol{u}_{i}$,
where $\boldsymbol{\beta}_{n1}$ is a $q_{n}\times 1$ vector corresponding to
significant covariate effects $\boldsymbol{X}_{i(1)}$ and
$\boldsymbol{\beta}_{n2}$ is a $(p_{n}-q_{n})\times 1$ vector of zeros, and
$g_{ik}(.)=0$ almost surely for $k=s_{n}+1,\ldots,r_{n}$, and $s_{n}$ is the
total number of nonzero function components. Therefore
$\\{\boldsymbol{\alpha}_{nk}\\}_{k=1}^{s_{n}}$ are nonzero vectors associated
to relevant functions $\\{g_{ik}(.)\\}_{k=1}^{s_{n}}$ and the
$\\{\boldsymbol{\alpha}_{nk}\\}_{k=s_{n}+1}^{r_{n}}$ are vectors of zero
corresponding to redundant functions $\\{g_{ik}(.)\\}_{k=s_{n}+1}^{r_{n}}$.
For technical convenience let
$\boldsymbol{\theta}_{n0}=(\boldsymbol{\theta}_{n01}^{\top},\boldsymbol{\theta}_{n02}^{\top})^{\top}$
be the true parameters, where
$\boldsymbol{\theta}_{n01}=(\boldsymbol{\beta}_{n01}^{\top},\boldsymbol{\alpha}_{n01}^{\top},\ldots,\boldsymbol{\alpha}_{n0s_{n}}^{\top})^{\top}$
is $(s=q_{n}+\sum_{k=1}^{s_{n}}h_{nk})$-dimensional vector of that the
elements are all nonzero and
$\boldsymbol{\theta}_{n02}=(\boldsymbol{\beta}_{n02}^{\top},\boldsymbol{\alpha}_{n0s_{n}+1}^{\top},\ldots,\boldsymbol{\alpha}_{n0r_{n}}^{\top})^{\top}=\boldsymbol{0}$
is a vector of dimension $(p_{n}-q_{n})+\sum_{k=s_{n}+1}^{r_{n}}h_{nk}$.
Consequently, estimated values and the design matrix is re-partitioned as
$\widehat{\boldsymbol{\theta}}_{n}=(\widehat{\boldsymbol{\theta}}_{n1}^{\top},\widehat{\boldsymbol{\theta}}_{n2}^{\top})^{\top}$,
and
$\boldsymbol{D}_{i}=\big{(}\boldsymbol{D}_{i(1)}^{\top},\boldsymbol{D}_{i(2)}^{\top}\big{)}^{\top}$
which
$\widehat{\boldsymbol{\theta}}_{n1}=(\widehat{\boldsymbol{\beta}}_{n1}^{\top},\widehat{\boldsymbol{\alpha}}_{n1}^{\top},\ldots,\widehat{\boldsymbol{\alpha}}_{ns_{n}}^{\top})^{\top}$,
$\widehat{\boldsymbol{\theta}}_{n2}=(\widehat{\boldsymbol{\beta}}_{n2}^{\top},\widehat{\boldsymbol{\alpha}}_{ns_{n}+1}^{\top},\ldots,\widehat{\boldsymbol{\alpha}}_{nr_{n}}^{\top})^{\top}$,
$\boldsymbol{D}_{i(1)}=\big{(}\boldsymbol{X}_{i(1)}^{\top},\boldsymbol{B}_{1}^{\top},\ldots,\boldsymbol{B}^{\top}_{s_{n}}\big{)}^{\top}$,
and
$\boldsymbol{D}_{i(2)}=\big{(}\boldsymbol{X}_{i(1)}^{\top},\boldsymbol{B}_{s_{n}+1}^{\top}\\\
,\ldots,\boldsymbol{B}_{r_{n}}^{\top}\big{)}^{\top}$. Throughout, we need some
regularity conditions and to save space we managed to put them in a separate
supplementary file (SF) and refer here as (A.1)-(A.7). Further, some lemmas
are also used that we put them in the SF. The following theorems establish the
selection and estimation consistency for both the functional regression
coefficients and regression parameters.

###### Theorem 1

Assume conditions (A.1)–(A.7) and that, the number of knots
$L_{n}=O_{p}(n^{1/(2r+1)})$ and the tuning parameter $\lambda_{n}\to 0$. Then
there exists a solution of
$\boldsymbol{U}_{n}(\boldsymbol{\theta}_{n})=o(a_{n})$, such that

$\displaystyle(i)$
$\displaystyle\|\widehat{\boldsymbol{\theta}}_{n}-\boldsymbol{\theta}_{n0}\|=O_{p}(\sqrt{p_{n}/n}),$
$\displaystyle(ii)$ $\displaystyle\max_{1\leq k\leq
p_{n}}\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{n_{i}}\big{(}\widehat{g}_{ik}(x_{ijk})-g_{0ik}(x_{ijk})\big{)}^{2}=O_{p}(n^{-2r/(2r+1)}).$

###### Theorem 2

Assume conditions (A.1)–(A.7), and that the tuning parameter $\lambda_{n}\to
0$, and $\lambda_{n}n^{r/(2r+1)}\to+\infty$. Then we have

$\displaystyle(i)$
$\displaystyle\mathbb{P}_{n}(\widehat{g}_{k}(.)=\boldsymbol{0})\to
1,\quad\forall~{}k=s_{n}+1,\ldots,p_{n},$ $\displaystyle(ii)$
$\displaystyle\mathbb{P}_{n}(\widehat{\boldsymbol{\beta}}_{n2}=\boldsymbol{0})\to
1,$ $\displaystyle(iii)$
$\displaystyle\boldsymbol{\xi}_{n}^{\top}\overline{\boldsymbol{M}}_{n}^{*^{-1/2}}(\boldsymbol{\beta}_{n0})\overline{\boldsymbol{H}}_{n}^{*}(\boldsymbol{\beta}_{n0})(\widehat{\boldsymbol{\beta}}_{n1}-\boldsymbol{\beta}_{n01})\overset{\mathcal{D}}{\to}\hbox{N}_{p_{n}}(0,1),$

where

$\displaystyle\overline{\boldsymbol{M}}_{n}^{*}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{u|y}\Big{[}\sum_{i=1}^{n}\boldsymbol{X}_{i}^{*^{\top}}\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\overline{\boldsymbol{R}}^{-1}\boldsymbol{R}_{0}\overline{\boldsymbol{R}}^{-1}\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\boldsymbol{X}_{i}^{*}\Big{]},$
$\displaystyle\overline{\boldsymbol{H}}_{n}^{*}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{u|y}\Big{[}\sum_{i=1}^{n}\boldsymbol{X}_{i}^{*^{\top}}\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\overline{\boldsymbol{R}}^{-1}\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\boldsymbol{X}_{i}^{*}\Big{]},$
$\boldsymbol{X}_{i}^{*}=(\boldsymbol{I}-\boldsymbol{P})\boldsymbol{X}_{i},\quad\boldsymbol{P}=\boldsymbol{B}(\boldsymbol{B}^{\top}\boldsymbol{\Omega}\boldsymbol{B})^{-1}\boldsymbol{B}^{\top}\boldsymbol{\Omega},\quad\boldsymbol{\Omega}=\hbox{diag}\\{\boldsymbol{\Omega}_{i}\\}$

and

$\boldsymbol{\Omega}_{i}=\mathbb{E}_{u|y}\Big{[}\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\overline{\boldsymbol{R}}^{-1}\boldsymbol{A}_{i}^{\frac{1}{2}}(\boldsymbol{\theta}_{n},\boldsymbol{u}_{i})\Big{]}.$

##  4 Numerical Studies

In this section, in a series of numerical expriments, we assess the
performance of our proposal.

###  4.1 Simulations

In this section, we conduct three simulation studies from the GAMM with both
continuous and binary outcomes for 100 simulated data sets. For evaluating the
estimation accuracy, we report the empirical mean square error (MSE), defined
as
$\sum_{r=1}^{100}\|\widehat{\boldsymbol{\beta}}_{n}^{r}-\boldsymbol{\beta}_{n0}\|/100$
where $\widehat{\boldsymbol{\beta}}_{n}^{r}$ is the estimator of
$\boldsymbol{\beta}_{n0}$ obtained using the $r$th generated data set. We
evaluate the total averaged integrated squared error (TAISE) to assess
estimation efficiency of the functional part. Let $\widehat{g}^{r}$ be the
estimator of a nonparametric function, $g$, in the $r$th replication and
$\\{x_{m}\\}_{m=1}^{n_{grid}}$ be the grid points where $\widehat{g}^{r}$ is
evaluated. We define

$\rm{AISE}(\widehat{g})=\frac{1}{100}\sum_{r=1}^{100}\frac{1}{n_{grid}}\sum_{m=1}^{n_{grid}}\\{\widehat{g}^{r}(x_{m})-g(x_{m})\\}^{2}$

and $\rm{TAISE}=\sum_{k=1}^{p_{n}}\rm{AISE}(\widehat{g}_{k})$. The performance
in variable selection is gauged by ’FZf’, number of false zero functional
predictors; ’FNf’, number of false nonzero functional predictors; ’FZs’,
number of false zero scalar covariates; ’FNs’, number of false nonzero scalar
covariates. To present a more comprehensive picture, we also use other
criteria for variable selection performance evaluation. ’U.fit’, ’C.fit’, and
’O.fit’ give the percentage of under fitting, correct fitting and over fitting
from 100 replications, respectively.

####  4.1.1 Example 1: Continuous response and moderate dimension of
covariates

In this example, the continuous responses $y_{ij}$ conditional on independent
random effect $u_{i}\sim\mathcal{N}(0,0.5)$ are generated from

$\displaystyle
y_{ij}|u_{i}=\sum_{k=1}^{r_{n}}x_{ijk}\beta_{k}+\sum_{k=1}^{p_{n}}g_{k}(x_{ijk})+u_{i}+\epsilon_{ij},\quad
i=1,\ldots,p_{n};\qquad j=1,\ldots,5,$ (9)

where $p_{n}=10$ and the number of sample size is $n=100,250$ or $500$. The
true regression coefficients are $\boldsymbol{\beta}=(-1,-1,2,0,\ldots,0)$ and
the additive functions are

$g_{1}(x_{ij1})=(2x_{ij1}-1)^{2},\quad g_{2}(x_{ij2})=8(x_{ij2}-0.5)^{3},\quad
g_{2}(x_{ij2})=\sin(2\pi x_{ij3})$

and $g_{k}(x_{ijk})=0$ for $k=4,\ldots,10$. Thus the last seven coefficients
and functional variables in this model are redundant and do not contribute to
the model. The covariates
$\boldsymbol{X}_{ij}=(x_{ij1},\ldots,x_{ijp_{n}})^{\top}$ are generated
independently from Uniform $([0,1]^{10})$. The error
$\boldsymbol{\epsilon}_{i}=(\epsilon_{i1},\ldots,\epsilon_{i5})^{\top}$
follows a multi-5-variate normal distribution with mean $0$, a common marginal
variance $\sigma^{2}=1$, and an AR-1 correlation structure with correlation
$\rho=0.7$.

We apply the PGEE with SCAD penalty and fit the linear splines ($d$ = 1) and
cubic splines ($d$ = 3). To illustrate the effect incorporating within-cluster
correlation on estimation efficiency, we compare the estimation efficiency of
using basis matrices from different working correlation structures:
exchangeable (EX), AR-1, and independent (IND). In addition, we compare the
PGEE approach with the GEE estimations of a full model (FULL) and an oracle
model (ORACLE). Here, the full model consists of all 10 variables and the
oracle model contains only the first three relevant variables.

Table 1 summarizes the estimation results of all procedures. It shows that the
estimation procedures with a correct AR-1 working correlation have the
smallest MSEs and TAISEs, and thus the estimators are more efficient than
their counterparts with IND working correlation that ignore within-cluster
correlation. Estimation based on a misspecified EX correlation structure will
lead to some efficiency loss compared with using the true AR-1 structure, but
it is still more efficient than assuming IND structure. In general, estimation
using cubic splines is more efficient than that using linear splines.
Furthermore, with the same type of working correlation, ORACLE has the
smallest MSEs and TAISEs compared with the SCAD and FULL model approaches.
This is not surprising, because ORACLE takes advantage of the true information
from the data generation and contains only three relevant variables. However,
the efficiency of the SCAD is closer to that of ORACLE, with much smaller MSEs
and TAISEs than those of the FULL model.

Table 2 provides variable selection results for the SCAD procedures. It shows
that the selection of functional and scalar predictors of the proposed method
is quite accurate and stable with a very small number of FNs. Also for the
AR-1, EX, or IND type of working correlation, the percentage of correct model
fitting goes to 1 quickly as the sample size increases. This confirms the
consistency theorem of variable selection provided in Section 3.2.
Furthermore, for the cubic spline approach, the correct modelfitting
percentage using the AR-1 correlation structure is noticeably higher than that
using the IND for all sample size.

Table 1: Example 1: Continuous response and moderate dimensional covariates.
Estimation results: the MSEs (TAISEs in parentheses) from 100 simulations.

---
n | method | linear spline | cubic spline
|  | AR-1 | EX | IND | AR-1 | EX | IND
100 | SCAD | 0.041(1.002) | 0.250(3.561) | 0.382(5.507) | 0.020(0.087) | 0.028(0.240) | 0.027(0.118)
ORACLE | 0.041(0.179) | 0.067(0.328) | 0.071(0.136) | 0.019(0.179) | 0.024(0.167) | 0.032(0.184)
FULL | 0.148(1.148) | 0.970(6.658) | 2.654(4.780) | 0.234(0.103) | 1.778(0.374) | 0.269(0.160)
250 | SCAD | 0.027(0.661) | 0.090(1.266) | 0.054(1.159) | 0.013(0.056) | 0.013(0.059) | 0.018(0.078)
ORACLE | 0.016(0.128) | 0.069(0.340) | 0.047(0.086) | 0.011(0.117) | 0.013(0.084) | 0.018(0.064)
FULL | 0.096(0.765) | 0.331(2.348) | 0.602(4.251) | 0.143(0.064) | 0.207(0.099) | 0.160(0.093)
500 | SCAD | 0.019(0.466) | 0.022(0.501) | 0.030(0.753) | 0.009(0.039) | 0.009(0.039) | 0.013(0.053)
ORACLE | 0.011(0.094) | 0.044(0.211) | 0.030(0.057) | 0.007(0.008) | 0.008(0.060) | 0.023(0.000)
FULL | 0.068(0.543) | 0.235(2.040) | 0.091(0.836) | 0.102(0.046) | 0.111(0.047) | 0.115(0.063)

Table 2: Example 1: Continuous response and moderate dimensional covariates.
Variable selection results from 100 simulations.

---
n | method |  | FZs | FNs | FZf | FNf | U.fit | C.fit | O.fit
100 | linear spline | AR-1 | 0.020 | 0.080 | 0.000 | 0.000 | 0.02 | 0.92 | 0.06
EX | 0.030 | 0.250 | 0.000 | 0.000 | 0.03 | 0.83 | 0.14
IND | 0.000 | 0.630 | 0.000 | 0.000 | 0.00 | 0.52 | 0.48
cubic spline | AR-1 | 0.000 | 0.080 | 0.000 | 1.000 | 0.00 | 0.92 | 0.08
EX | 0.070 | 0.020 | 0.000 | 1.000 | 0.07 | 0.92 | 0.01
IND | 0.060 | 0.260 | 0.000 | 1.000 | 0.06 | 0.70 | 0.24
250 | linear spline | AR-1 | 0.020 | 0.000 | 0.000 | 0.000 | 0.02 | 0.98 | 0.00
EX | 0.080 | 0.050 | 0.000 | 0.000 | 0.08 | 0.89 | 0.03
IND | 0.270 | 0.000 | 0.000 | 0.000 | 0.27 | 0.73 | 0.00
cubic spline | AR-1 | 0.000 | 0.000 | 0.000 | 1.000 | 0.00 | 1.00 | 0.00
EX | 0.040 | 0.000 | 0.000 | 1.000 | 0.04 | 0.96 | 0.00
IND | 0.180 | 0.020 | 0.000 | 1.000 | 0.18 | 0.80 | 0.02
500 | linear spline | AR-1 | 0.000 | 0.000 | 0.000 | 0.000 | 0.00 | 1.00 | 0.00
EX | 0.040 | 0.430 | 0.000 | 0.000 | 0.04 | 0.78 | 0.18
IND | 0.170 | 0.030 | 0.000 | 0.000 | 0.17 | 0.80 | 0.03
cubic spline | AR-1 | 0.000 | 0.020 | 0.000 | 1.000 | 0.00 | 0.98 | 0.02
EX | 0.000 | 0.000 | 0.000 | 1.000 | 0.00 | 1.00 | 0.00
IND | 0.130 | 0.000 | 0.000 | 1.000 | 0.13 | 0.87 | 0.00

####  4.1.2 Example 2: Continuous response with high-dimensional covariates

To assess our method in more challenging cases for high-dimensional data, we
consider a model with the dimension of functional components $p_{n}=100$ in
(9). However, only the first three variables are relevant and take the same
functional forms as in Example 1. We consider the model (9) with moderate
number of $n=200$ and $n_{i}=5$, with errors
$\\{\boldsymbol{\epsilon}_{i}=(\epsilon_{i1},\ldots,\epsilon_{i5})^{\top}\\}_{i=1}^{200}$
generated independently from a multivariate normal distribution with mean $0$,
and an AR-1 correlation structure with $\mathop{\rm
corr}\nolimits(\epsilon_{ij},\epsilon_{ij^{{}^{\prime}}})=0.7^{|j-j^{{}^{\prime}}|}$
for $1\leq j,j^{{}^{\prime}}\leq 5$. We apply the GEE to estimate the full and
oracle models with working independent, exchangeable, or AR-1 working
correlation. For variable selection, we consider the PGEE with SCAD penalty
with basis matrices from IND, EX, or AR-1 working correlations. Table 3
reports MSEs and TAISEs for FULL, ORACLE, and SCAD and variable selection
results on correct, overfit, and underfit percentages of the SCAD approach for
three IND, EX and AR-1 working correlations. Table 3 clearly indicates that
the improvement from incorporating within-cluster correlation is very
significant. In particular, the estimation procedures with a correctly
specified AR-1 structure always give smaller MSEs and TAISEs than those with a
misspecified EX or IND working correlation. Also, the efficiency gained by
incorporating correlation could be increased in the cubic splines approach.
For variable selection, the SCAD in both cases of linear and cubic, with an
AR-1 working correlation also performs noticeably better than the one with EX
or IND working correlation. Furthermore, Table 3 also shows that the SCAD
procedure dramatically improves the estimation accuracy for this high-
dimensional case, with smaller MSEs and TAISEs from the FULL model.

Table 3: Example 2: Continuous response and high-dimensional covariates. The
MSEs (TAISEs in parentheses), percentages of correct fitting (C.fit),
underfitting (U.fit), and overfitting (O.fit) from 100 replications.

---
|  | SCAD | ORACLE | FULL | FZs | FNs | FZf | FNf | U.fit | C.fit | O.fit
linear spline | AR-1 | 0.018(0.073) | 0.027(0.659) | 1.662(0.576) | 0.02 | 0.00 | 0.00 | 1.00 | 0.02 | 0.98 | 0.00
EX | 0.031(0.113) | 0.056(1.003) | 1.056(0.837) | 0.01 | 2.22 | 0.00 | 0.00 | 0.01 | 0.94 | 0.05
IND | 0.029(0.197) | 0.078(1.002) | 3.206(4.991) | 0.18 | 0.02 | 0.00 | 0.00 | 0.17 | 0.81 | 0.02
cubic spline | AR-1 | 0.014(0.062) | 0.013(0.056) | 2.273(0.081) | 0.10 | 0.00 | 0.00 | 0.00 | 0.07 | 0.93 | 0.00
EX | 0.016(0.088) | 0.014(0.152) | 2.571(0.142) | 0.15 | 0.00 | 0.00 | 1.00 | 0.12 | 0.88 | 0.00
IND | 0.011(0.147) | 0.018(0.076) | 3.017(4.617) | 0.27 | 0.35 | 0.00 | 1.00 | 0.26 | 0.68 | 0.06

####  4.1.3 Example 3: Binary response

To assess the performance of our method for binary outcomes, we generate a
moderate random sample of 200 subjects in each simulation run. Within each
subject, binary responses $\\{y_{ij}\\}_{j=1}^{200}$ are generated from a
marginal logit model

$\mathop{\rm
logit}\nolimits\mathbb{P}(y_{ij}=1|\boldsymbol{X}_{ij}=\boldsymbol{X}_{ij})=\sum_{k=1}^{10}x_{ijk}\beta_{k}+\sum_{k=1}^{10}g_{k}(x_{ijk})+u_{i}$

where $g_{1}(x_{ij1})=[\exp(x+1)-\big{(}\exp(2)-\exp(1)\big{)}/16]$,
$g_{2}(x_{ij2})=\cos(2\pi x_{ij2})/4$ $g_{3}(x_{ij3})=x_{ij3}(1-x_{ij3})-1/6$
$g_{4}(x_{ij4})=2(x_{ij4}-0.5)^{3}$ and the remaining $6$ covariates are null
variables with $g_{4}(x_{ij4})=0$ for $k=5,\ldots,10$. The covariates
$\boldsymbol{X}_{ij}=(x_{ij1},\ldots,x_{ijp_{n}})^{\top}$ are generated
independently from Uniform $([0,1]^{10})$. We generate correlated binary
responses with AR-1 correlation structure with a correlation parameter of
$0.5$. We applied the PGEE with a SCAD penalty for the variable selection, and
the GEE for estimation of the full and oracle models. To illustrate how
different working correlations could affect our estimation and variable
selection results, we consider EX and IND structures, in addition to the true
AR-1 correlation structure. Table 4 gives the TAISEs for the SCAD, ORACLE, and
FULL models with three different working correlations. Similar to the previous
continuous simulation studies, MSEs and TAISEs calculated based on the both
linear and cubic SCAD approach are also shown to be close to the MSEs and
TAISEs from ORACLE, and much smaller than those from the FULL model.

For all three aformentioned study, the amounts of FZs, FNs, FZf, FNf, U.fit,
C.fit and O.fit shows that the SCAD approach can gain significant estimation
accuracy by effectively removing the zero component variables. Overall, the
SCAD procedures work reasonably well, and the SCAD with EX and AR-1 working
correlation structures provides better variable selection results than the
SCAD with IND working structure. In addition, the SCAD approach with AR-1
correlation performs the best in selecting non-zero component variables.

Table 4: Example 3: Binary response. The MSEs (TAISEs in parentheses),
percentages of correct fitting (C.fit), underfitting (U.fit), and overfitting
(O.fit) from 100 replications.

---
|  | SCAD | ORACLE | FULL | FZs | FNs | FZf | FNf | U.fit | C.fit | O.fit
linear spline | AR-1 | 0.443(0.249) | 0.209(0.000) | 0.456(0.972) | 0.00 | 0.00 | 0.00 | 0.00 | 0.000 | 1.00 | 0.00
EX | 0.451(0.262) | 0.207(0.000) | 0.456(0.988) | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | 0.00
IND | 0.447(0.236) | 0.209(0.000) | 0.455(0.980) | 0.01 | 0.00 | 0.00 | 0.00 | 0.01 | 0.99 | 0.00
cubic spline | AR-1 | 0.474(0.000) | 0.257(0.000) | 1.899(1.759) | 0.14 | 0.00 | 0.00 | 2.00 | 0.14 | 0.86 | 0.00
EX | 0.478(0.000) | 0.257(0.000) | 1.903(1.761) | 0.15 | 0.00 | 0.00 | 2.00 | 0.15 | 0.85 | 0.00
IND | 0.477(0.000) | 2.579(0.000) | 1.882(1.746) | 0.15 | 0.00 | 0.00 | 2.00 | 0.15 | 0.85 | 0.00

###  4.2 Benchmark analysis

In this subsection, we analyse the data from the Multi-Center AIDS Cohort
study. The data set contains the human immunodeficiency virus, HIV, status of
283 homesexual men infected with HIV during the follow-up period between 1984
to 1991. All individual were scheduled to have their measurment made during
semi annual visit. Here $t_{ij},i=1,\ldots,n,j=1,\ldots,n_{i}$ denote the time
length in years between seroconversion and the $j$th measurment of the $i$th
individual after the infection. Xue and Zhu (2007) analyzed the data set using
partial linear models. The primary interest was to describe the trend of the
mean CD4 percentage depletion over time and to evaluate the effect of
cigarette smoking, pre-HIV infection CD4 percentage, and age at infection on
the mean CD4 cell percentage after the infection. We take four covariates for
this study and construct an additive model by including all these covariates
in both parametric and functional parts. We apply the PGEE with SCAD penalty
and fit the linear splines ($d$ = 1) and cubic splines ($d$ = 3). In addition,
we compare the PGEE approach with the GEE estimations of a full model. We use
the standard errors (SE) were all calculated using the leave-one-out cross
validation method. To best identify a model supported by the data, mean
squared estimation errors (MSEE) and mean squared prediction errors (MSPE)
were calculated for each method.

Table 5 reports the results. Two variable selection methods resulted in models
with smaller SE, MSEE, and better predictions than the full model, which
indicates that redundant variables exist in data sets. Cubic spline gave
parsimonious models with the best prediction performance. Linear spline select
the all nonparametric functions as zero, but cubic spline and FULL models
select time variable is nonparametric. As shown in figure 1, we see that mean
baseline CD4 perecentage of the population (blue line) and all subjects (gray
lines) depreters rather quickly at the begining of HIV infection, but the CD4
perecentage appears to be incresing two years after the infection. Maybe
taking medication has caused such behavior.

Table 5: CD4 data example: estimation, prediction and variable selection
results.

---
| linear spline | cubic spline | FULL
variabels | estimete | SD | estimete | SD | estimete | SD
time | 0.000 | 0.000 | -6.247 | $1.4\times 10^{-4}$ | -0.448 | 1.069
pre-CD4 | -0.342 | 0.114 | 0.175 | $3.5\times 10^{-4}$ | -0.457 | 0.095
age | 0.556 | 0.126 | 1.576 | $5.3\times 10^{-4}$ | 0.359 | 0.113
smoke | 0.000 | 0.000 | -0.953 | $4.2\times 10^{-6}$ | -0.515 | 1.512
MSEE | 0.240 | 0.049 | 13.623
MSPE | 10.101 | 9.422 | 11.534

Figure 1: The estimated component function ($\widehat{f}(t)$) from the SCAD
procedure using the cubic spline.

###  4.3 Yeast Cell-Cycle Gene Expression data analysis

A yeast cell-cycle gene expression data collected in the CDC15 experiment of
Spellman et al. (1998) where genome-wide mRNA levels of 6178 yeast ORFs
(abbreviation for open reading frames, which are DNA sequences that can
determine which amino acids will be encoded by a gene) at 7 minute intervals
for 119 minutes, which covers two cell-cycle periods for a total of 18 time
points, measured. The cell cycle is a tightly regulated life process where
cells grow, replicate their DNA, segregate their chromosomes, and divide into
as many daughter cells as the environment allows. The cell-cycle process is
commonly divided into M/G1-G1-S-G2-M stages. Refer to Wang et al. (2012), for
more detailed description of this data set. Transcription factors (TFs) have
been observed to play critical roles in gene expression regulation. A TF
(sometimes called a sequence-specific DNA-binding factor) is a protein that
binds to specific DNA sequences, thereby controlling the flow (or
transcription) of genetic information from DNA to mRNA. To better understand
the phenomenon underlying cell-cycle process, it is important to identify TFs
that regulate the gene expression levels of cell cycle-regulated genes. It is
not clear where these TFs regulate all cell cycle genes, however. We applied
our methods to identify the key TFs. The dataset that we use present a subset
of 283 cell-cycled-regularized genes observed over 4 time points at G1 stage.
The response variable $Y_{ij}$ is the log-transformed gene expression level of
gene $i$ measured at time point $j$, for $i=1,\ldots,283$. We use the
following GAMM model

$\displaystyle
y_{ij}=\sum_{k=1}^{96}x_{ijk}+\sum_{k=1}^{96}g_{ik}(x_{ijk})+u_{i},$

where the covariates $x_{ijk}$ , $k=1,\ldots,96$, is the matching score of the
binding probability of the $k$th TF on the promoter region of the $i$th gene.
The binding probability is computed using a mixture modeling approach based on
data from a ChIP binding experiment; see Wang et al. (2007) for details.
Covariates $x_{ijk}$ is standardized to have mean zero, and $u_{i}$ is the
random intercept. Our goal is to identify the TFs that might be related to the
expression patterns of these 283 cellcycle–regulated genes. Therefore we apply
a penalization procedure by the proposal P-GAMM and fit the linear splines
($d$ = 1) and cubic splines ($d$ = 3), and also use the full model for
comparing. Table 6 reports the estimation and prediction performance when
linear and cubic splines and full model are adopted. For stage G1, MBP1, SWI4,
and SWI6 are three TFs that have been proved important in the aforementioned
biological experiments and our analysis reveals that they have been selected
by the penalization method with linear and cubic splines. Among this important
TFs, linear spline select the all nonparametric effects as zero, but cubic
spline identify that MBP1 and SWI4 have nonparametric effects too. Our
proposal in both cases of linear and cubic give smaller MSEE and MSPE then the
FULL model, where show that the penalized GAMM improves the estimation and
prediction performance. Prediction accuracy are similar in both linear and
cubic splines, but linear spline perform more efficient in term of MSEE.

Table 6: Yeast Cell-Cycle Gene Expression data: estimation and prediction
performance.

---
| linear spline | cubic spline | FULL
MSEE | 0.002 | 0.176 | 4.141
MSPE | 0.351 | 0.339 | 0.463

##  5 Concluding Remarks

In general, when the number of covariates is large, identifying the exact
underlying model is a challenging task, in particular when some of the nonzero
signals are relatively weak. Here, we consider the GAMM when longitudinal
responses collected in a high dimensional regim. In our setting, there is no
specified likelihood function for the GAM, because the outcomes could be
nonnormal and discrete, which makes estimation and model selection very
challenging problems. Also, we found that it may be quite computationally
intensive in high-dimensional variable selection settings, because the
dimension of the parameters in the nonparametric forms increases significantly
compared with parametric model selection settings. We approximated the
additive functional components based on polynomial B-spline smoothing and
propose doubel penalized method where variable selection and estimation of the
both parametric and nonparametric components obtained, and were able to select
the functional components groupwise. The procedure involved the specification
of the posterior distribution of the random effects, which cannot be evaluated
in a closed form. We used a Metropolis algorithm which does not require the
specification of the posterior distribution. To implement the procedure in
practice, a computationally flexible iterative algorithm developed. We
established an asymptotic property with consistency for the fixed effect
covariates and nonparametric components, which achieves the optimal rate of
convergence. In addition, the proposed model selection strategy was able to
select the correct GAMM consistently. That is, with probability approaching to
1, the estimators for the zero function components converge to 0 almost
surely. We illustrated our method using numerical studies with both continuous
and binary responses, along with real data applications. Results demonstrated
that the proposal works well, and can correctly select the nonzero fixed
effect covariates and functional components with probability tending to one as
the sample size goes to infinity.

## supplement

The regularity conditions (A.1)-(A.7), proofs of the main results, and some
instrumental lemmas are provided in a separate supplementary file.

## References

  * Berhane and Tibshirani (1998) Berhane, K., and Tibshirani, R. J. (1998). Generalized additive models for longitudinal data. Canadian Journal Statistics. 26, 517–535.
  * Breslow and Clayton (1993) Breslow, N. E., and Clayton, D. G. (1993). Approximate inference in generalized linear mixed models. J. Amer. Statist. Assoc. 88 9–25.
  * Carroll et al. (2009) Carroll, R. J., Maity, A., Mammen, E., and Yu, K. (2009). Nonparametric additive regression for repeatedly measured data. Biometrika, 96 383–398.
  * Chu et al. (2020) Chu, W., Li, R., Liu, J., and Reimherr, Matthew. (2020). Feature selection for generalized varying coefficient mixed-effect models with application to obesity GWAS. Ann. Appl. Stat. 14 (1) 276–298.
  * Chu et al. (2016) Chu, W., Li, R., and Reimherr, M. (2016). Feature screening for time-varying coefficientmodelswith ultrahigh dimensional longitudinal data. Ann. Appl. Statist. 10 596–617.
  * Dziak (2006) Dziak, J. J. (2006). Penalized quadratic inference functions for variable selection in longitudinal research. Ph.D Thesis, the Pennsylvania State University. (https://etda.libraries.psu.edu/paper/7084/)
  * Fan and Li (2001) Fan, J. Q., and Li, R. (2001). Variable selection via nonconcave penalized likelihood and its oracle properties. J. Amer. Statist. Assoc. 96 1348–1360.
  * Fan and Li (2004) Fan, J. Q., and Li, R. (2004). New estimation and model selection procedures for semiparametric modeling in longitudinal data analysis. J. Amer. Statist. Assoc. 99 710–723.
  * Fan and Li (2012) Fan, Y., and Li, R. (2012). Variable selection in linear mixed effects models. Ann. Statist. 40(4) 2043–2068.
  * Fang et al. (2020) Fang, X. E., Ning, Y., and Li. R (2020). Test of significance for high-dimensional longitudinal data. Ann. Statist. 48(5) 2622–2645
  * Fu (2003) Fu, W. J. (2003). Penalized estimating equations. Biometrics 59 126–132.
  * Groll and Tutz (2012) Groll, A., and Tutz. G. (2012). Regularization for generalized additive mixed models by likelihood-based boosting. Methods Inf Med. 51 168–177.
  * Li and Liang (2008) Li, R., and Liang, Hua. (2008). Variable selection in semiparametric regression modeling. Ann. Statist. 36(1) 261–286.
  * Lin and Zhang (1999) Lin, x., and Zhang, D. (1999). Inference in generalized additive mixed models by using smoothing spline. J. Royal Statist. Soc. Ser. B 61 381–400.
  * Ma et al. (2013) Ma, S., Song, Q., and Wang, L. (2013). Simultaneous variable selection and estimation in semiparametric modeling of longitudinal/clustered data. Bernoulli 19 252–274.
  * Ni et al. (2010) Ni, X., Zhang, D., and Zhang, H. H. (2010). Variable selection for semiparametric mixed models in longitudinal studies. Biometrics 66 79–88.
  * Spellman et al. (1998) Spellman, P. T., Sherlock, G., Zhang, M. Q., Iyer, V. R., Anders, K., Eisen, M. B., Brown, P. O., Botstein, D., and Futcher, B. (1998). Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization. Molecular Biology of Cell, 9, 3273–3297.
  * Tanner (1993) Tanner, M. A. (1993). Tools for Statistical Inference: Methods for the Exploration of Posterior Distributions and Likelihood Functions (2nd ed.), New York: Springer-Verlag.
  * Taavoni and Arashi (2021) Taavoni, M., and Arashi, M. (2021). High-dimensional generalized semiparametric model for longitudinal data. Statistics A J. Theo. App. Statist. 55(4), 831–850.
  * Wang et al. (2007) Wang, L., Chen, G., and Li, H. (2007). Group SCAD regression analysis for microarray time course gene expression data. Bioinformatics 23 1486–1494.
  * Wang et al. (2008) Wang, L., Li, H., and Huang, J. Z. (2008). Variable selection in nonparametric varying-coefficient models for analysis of repeated measurements. J. Amer. Statist. Assoc. 103 1556–1569.
  * Wang et al. (2012) Wang, L., Zhou, J., and Qu, A. (2012). Penalized generalized estimating equations for high-dimensional longitudinal data analysis. Biometrics 68 353–360.
  * Wu and Xue (2014) Wu, J. and Xue, L. (2014). Model detection for additive models with longitudinal data. Open J. Statist. 4 868–878.
  * Xue (2009) Xue., L. (2009). Consistent variable selection in additive models. Statistica Sinica. 19 1281–1296.
  * Xue and Liang (2009) Xue, L., and Liang, H. (2009). Polynomial spline estimation for A generalized additive coefficient model. Scand. Statist. Theory Appl. 37 26–46.
  * Xue et al. (2010) Xue, L., Qu, A., and Zhou, J. (2010). Consistent model selection for marginal generalized additive model for correlated data. J. Amer. Statist. Assoc. 105 1518–1530.
  * Xue and Zhu (2007) Xue, L. G., and Zhu, L. X. (2007). Empirical likelihood for a varying coefficient model with longitudinal data. J. Amer. Statist. Assoc. 102 642–654.
  * You and Zhou (2007) You, J. and Zhou, H. (2007). Two-stage efficient estimation of longitudinal nonparametric additive models. Statistics and Probability Letters, 77 1666–1675.
  * Zeger and Karim (1991) Zeger, S. L., and Karim, M. R. (1991). Generalized linear models with random effects: A Gibbs sampling approuch. J. Amer. Statist. Assoc. 86 79–86.


counter-example.

Both claims have been machine-checked. For the $6$-element counter-example,
$R$ isn’t anti-symmetric due to $aRb\land bRa$, not semi-connex since $a$ and
$c$ are incomparable, and not dense since $dRe$ has no intermediate element;
left and right seriality has been achieved by making the elements at the start
and at the end of each chain reflexive.

Figure 34: Relation graph in Exm. 89

###### Example 89

[1]On a set $X$ of $6$ elements, a relation $R$ can satisfy semi-order
property 1 and be semi-connex, but neither quasi-transitive nor dense nor left
nor right serial (269). An example is shown in Fig. 34. It is not quasi-
transitive, since $eRc\land\lnot cRe$ and $cRf\land\lnot fRc$, but $fRe$; it
is not dense, since $bRa$ has no intermediate element; it is not left and
right serial, since $d$ and $a$ has no predecessor and successor,
respectively. On a set of $\leqslant 5$ elements, no relation with these
properties exists.

Figure 35: Relation graph in Exm. 90

###### Example 90

[1]On a set $X$ of $\leqslant 5$ elements, an irreflexive, semi-connex, and
dense relation $R$ must be quasi-transitive, or left or right serial (270). On
a larger set, this is no longer true; Fig. 35 shows a counter-example for the
$6$-element set $X=\\{a,b,c,d,e,f\\}$,

Both claims have been machine-checked. In Fig. 35, the cycle $a$, $b$, $c$
(shown in blue) violates quasi-transitivity, the minimal and maximal element
$d$ and $e$ (green and red) violates left and right seriality, respectively.
Element $f$ (grey) is related to each of $a,b,c$ in both directions, e.g.
$aRf\land fRa$ holds, thus helping to establish density.

Figure 36: Relation graph in Exm. 91

###### Example 91

[1]On a set $X$ of $\leqslant 5$ elements, an transitive and incomparability-
transitive relation $R$ must be anti-symmetric, semi-connex, dense, or left or
right serial (271). On a larger set, this is no longer true; Fig. 36 shows a
counter-example for the $6$-element set $X=\\{a,b,c,d,e,f\\}$, It is not anti-
symmetric, since $aRb$ and $bRa$; it is not semi-connex, since $d,f$ are
incomparable; it is not dense, since $cRf$ has no intermediate element; it is
not left and right serial, since $e$ and $f$ has no predecessor and successor,
respectively.

Figure 37: Relation graph in Exm. 92

###### Example 92

[1]On a set $X$ of $6$ elements, a relation $R$ can satisfy semi-order
properties 1 and 2 and be left serial but neither right serial nor
incomparability-transitive nor quasi-transitive nor dense (272). An example is
shown in Fig. 37. It is not right serial, since $b$ has no successor; it is
not incomparability-transitive, since $b,a$ and $b,d$ are incomparable, but
$dRa$; it is not quasi-transitive, since $dRf\land\lnot fRd$ and
$fRc\land\lnot cRf$, but not $dRc$; it is not dense, since $cRb$ has no
intermediate element. By reverting the arrow directions, a counter-example for
the dual 273 is obtained. On a set of $\leqslant 5$ elements, no relation with
either property combination exists.

Figure 38: Relation graph in Exm. 93

###### Example 93

[1]On a set $X$ of $6$ elements, a relation $R$ can satisfy semi-order
properties 1 and 2 and be quasi-transitive and left and right serial, but
neither transitive nor incomparability-transitive nor dense (274). An example
is shown in Fig. 38. It is not transitive, since $cRf\land fRc$ but not $cRc$;
it is not incomparability-transitive, since $a,b$ and $a,d$ are incomparable,
but $dRb$; it is not dense, since $dRb$ has no intermediate element. On a set
of $\leqslant 5$ elements, no relation with these properties exists.

##  6 Implementation issues

In this section, we comment on some program details. The source code is
available in the ancillary file nonprominentProperties.c.

###  6.1 Improved relation enumeration

$\begin{array}[]{l|cccc|l @{\hspace*{2cm}}
l|cccc|l}\lx@intercol\hfil\lx@intercol&a&b&c&\lx@intercol\hfil
d\hfil\lx@intercol&\hfil\hskip
56.9055pt&\lx@intercol\hfil\lx@intercol&b&c&d&\lx@intercol\hfil
a\hfil\lx@intercol&\\\ \cline{2-5}\cr\cline{8-11}\cr a&1&1&1&1&\langle
3,1\rangle\hfil\hskip 56.9055pt&b&1&1&.&.&\langle 1,1\rangle\\\
b&.&1&1&.&\langle 1,1\rangle\hfil\hskip 56.9055pt&c&.&.&1&1&\langle
2,0\rangle\\\ c&1&.&.&1&\langle 2,0\rangle\hfil\hskip
56.9055pt&d&.&1&.&1&\langle 2,0\rangle\\\ d&1&.&1&.&\langle
2,0\rangle\hfil\hskip 56.9055pt&a&1&1&1&1&\langle 3,1\rangle\\\
\cline{2-5}\cr\cline{8-11}\cr\end{array}$

Figure 39: Example relation (left) and its normal form (right)

We encoded a binary relation by a square array of bools. Originally, we
enumerated all possible assignments of such an array. However, if $R$ is a
relation on a finite set $X$ of $n$ elements, and $\pi:X\rightarrow X$ is a
permutation, then $R^{\prime}$ defined by
$xR^{\prime}y:\Leftrightarrow\pi_{x}R\pi_{y}$ shares all properties from Def.
1 with $R$; that is, $R^{\prime}$ is reflexive iff $R$ is, etc.

In order to speed up relation enumeration, we therefore defined a normal form
for a square array as follows: To the $i$th array row, assign the pair
$\langle c_{i},d_{i}\rangle$, where $c_{i}$ is the number of true values in
columns $j\neq i$, and $d_{i}$ is the value of the $i$th column. An array is
in normal form if ${{\langle
c_{1},d_{1}\rangle}\leqslant\ldots\leqslant{\langle c_{n},d_{n}\rangle}}$
holds, where “$\leqslant$” denotes the lexicographic order.

Every array can be converted into normal form by applying a row-sorting
permutation simultaneous to rows and columns. Note that the pairs are designed
to be invariant under simultaneous row and column permutation.

As an example, the left half of Fig. 39 shows a relation on the set
$X=\\{a,b,c,d\\}$, and the pairs assigned to each row. For readability, we
denoted the value false and true by “.” and “$1$”, respectively. The right
half shows the corresponding normal form, obtained by sorting the rows by
ascending associated pairs, and permuting the columns in the same way.

    const bool allRows[72] = {
      /*  0: ***** row position 0 ********************* */
      /*  0: group <0,0> */   0,0,0,
      /*  3: group <0,1> */   1,0,0,
      /*  6: group <1,0> */   0,1,0, 0,0,1,
      /* 12: group <1,1> */   1,1,0, 1,0,1,
      /* 18: group <2,0> */   0,1,1,
      /* 21: group <2,1> */   1,1,1,
      /* 24: ***** row position 1 ********************* */
      /* 24: group <0,0> */   0,0,0,
      /* 27: group <0,1> */   0,1,0,
      /* 30: group <1,0> */   1,0,0, 0,0,1,
      /* 36: group <1,1> */   1,1,0, 0,1,1,
      /* 42: group <2,0> */   1,0,1,
      /* 45: group <2,1> */   1,1,1,
      /* 48: ***** row position 2 ********************* */
      /* 48: group <0,0> */   0,0,0,
      /* 51: group <0,1> */   0,0,1,
      /* 54: group <1,0> */   1,0,0, 0,1,0,
      /* 60: group <1,1> */   1,0,1, 0,1,1,
      /* 66: group <2,0> */   1,1,0,
      /* 69: group <2,1> */   1,1,1,
      /* 72: ***** end ******************************** */
    };
    const int gpS[3*2+1] = {  0,  3,  6, 12, 18, 21, 24 };
    int gp[card];
    int rw[card];

    void initFromRow(bool R[card][card]) {
      for (int r=0; r<card; ++r)
        for (int c=0; c<card; ++c)
          R[r][c] = allRows[ r*24 + rw[r] + cc ];
    }

    void check03(void) {
      bool R[card][card];
      for (gp[0]=0; gp[0]<card*2; ++gp[0])
        for (gp[1]=0; gp[1]<=gp[0]; ++gp[1])
          for (gp[2]=0; gp[2]<=gp[1]; ++gp[2])
            for (rw[0]=gpS[gp[0]]; rw[0]<gpS[gp[0]+1]; rw[0]+=card)
              for (rw[1]=gpS[gp[1]]; rw[1]<gpS[gp[1]+1]; rw[1]+=card)
                for (rw[2]=gpS[gp[2]]; rw[2]<gpS[gp[2]+1]; rw[2]+=card) {
                  initFromRow(R);
                  if (isRgEucl(R) && ! isTrans(R))
                    printRel(R);
                }
    }

Figure 40: Improved code to search for right Euclidean non-transitive
relations

It is sufficient to consider property combinations only for relation arrays
that are in normal form, and this is what our improved algorithm does. Figure
40 shows the improved code to search for right Euclidean non-transitive
relations over a $3$-element universe.111111 The unimproved code was shown in
Fig. 2 The list allRows contains all possible rows for a $3\times 3$ array
representing a relation. The rows are grouped by assigned pair; we have one
copy for each row position, to account for the varying column position of the
diagonal element. For example, in horizontal position $0$, the rows 1,1,0 and
1,0,1 are assigned the pair $\langle 1,1\rangle$; they are found at starting
index 12 of allRows. The corresponding row values for horizontal position $1$
are 1,1,0 and 0,1,1; they start at index 36. Procedure check03 iterates in the
loops on gp[0], gp[1], and gp[2] over all combinations of groups that lead to
a normal form, and in the loops on rw[0], rw[1], and rw[2], over all
combinations of rows from the current groups. The lists allRows and gpS were
precomputed by another program; its source code is available in the ancillary
file genTables.c.

Figure 41 shows, for set cardinalities $1$ through $7$, the number of all
array assignments (column “Unpruned”), the number of arrays in normal form
(column “Pruned”), and the quotient of both numbers, indicating the speed-up
factor. Figure 42 shows, for each property, the number of satisfying relations
found with the old (column “Old”) and with the improved (column “Pruned”)
enumeration method.

Card | Unpruned | Pruned | Ratio
---|---|---|---
1 | 2 | 2 | 1.0
2 | 16 | 10 | 1.6
3 | 512 | 140 | 3.6
4 | 65 536 | 6 170 | 10.6
5 | 33 554 432 | 907 452 | 36.9
6 | 68 719 476 736 | 460 631 444 | 149.1
7 | 562 949 953 421 312 | 827 507 617 792 | 680.2
Figure 41: Number of relations vs. carrier set cardinality

###  6.2 Quine-McCluskey implementation

The procedure computeLaws iterates over all relations, determining for each
the set121212 encoded as bit vector, see Fig. 42 of its properties, and
counting the number of occurrences of each such vector. After that, it calls
the Quine-McCluskey implementation qmc to compute all prime implicants of the
non-occurring vectors. The latter procedure performs a top-down breadth-first
search on the search graph.

Count | Count | Name | Def. | Encoding
---|---|---|---|---
(Old) | (Pruned) |  |  |
1 | 1 | Empty |  | `0x000001`
1 | 1 | Univ |  | `0x000002`
32 | 6 | CoRefl |  1.3 | `0x000004`
3 163 | 166 | LfEucl |  1.16 | `0x000008`
3 163 | 186 | RgEucl |  1.15 | `0x000010`
7 776 | 440 | LfUnique |  1.23 | `0x000020`
7 776 | 1 818 | RgUnique |  1.24 | `0x000040`
32 768 | 1 012 | Sym |  1.7 | `0x000080`
47 462 | 4 841 | AntiTrans |  1.13 | `0x000100`
59 049 | 3 870 | ASym |  1.8 | `0x000200`
59 049 | 3 870 | Connex |  1.11 | `0x000400`
154 303 | 3 207 | Trans |  1.12 | `0x000800`
467 750 | 11 103 | SemiOrd1 |  1.17 | `0x001000`
1 048 576 | 70 436 | Irrefl |  1.2 | `0x002000`
1 048 576 | 70 436 | Refl |  1.1 | `0x004000`
1 069 742 | 71 198 | QuasiRefl |  1.6 | `0x008000`
1 889 568 | 50 480 | AntiSym |  1.9 | `0x010000`
1 889 568 | 50 480 | SemiConnex |  1.10 | `0x020000`
3 756 619 | 113 142 | IncTrans |  1.22 | `0x040000`
4 498 393 | 144 128 | SemiOrd2 |  1.18 | `0x080000`
5 531 648 | 131 994 | QuasiTrans |  1.14 | `0x100000`
15 339 497 | 425 854 | Dense |  1.21 | `0x200000`
28 629 151 | 764 962 | LfSerial |  1.20 | `0x400000`
28 629 151 | 817 185 | RgSerial |  1.19 | `0x800000`
Figure 42: Number of relations on a $5$-element set

An example graph, showing all possible prime implicants for a Boolean function
of 3 variables is given in Fig. 43. At each node of the search graph, the
corresponding conjunction is checked by the procedure qmcRect: if no
combination in its covered rectangle131313 This terminology is inspired by the
Karnaugh diagram method; in Fig. 43, the rectangle covered by a node
corresponds to the set of all leaves below it.  is “off” and at least one is
“on”,141414 Since we are interested in _non_ -occurring vectors, “on”
corresponds to an occurrence count of zero, and “off” to a count $>0$. We
encode “don’t care” by a count of $-1$.  then it is actually a prime
implicant. In that case, we output its description using qmcPrint,151515 In
particular, we don’t perform the usual search of a _minimal_ set of prime
implicants covering all “on” vectors.  and set all vectors in its covered
rectangle to don’t care.

Note that we can’t perform a depth-first search: for example, if $a$ isn’t a
prime implicant, we can’t check its child $ab$ next, since it could satisfy
the above primeness criterion, but nevertheless be covered by a simpler prime
implicant, such as $b$.

${\it
true}$$a$$\overline{a}$$b$$\overline{b}$$c$$\overline{c}$$ab$$a\overline{b}$$ac$$a\overline{c}$$\overline{a}b$$\overline{a}\overline{b}$$\overline{a}c$$\overline{a}\overline{c}$$bc$$b\overline{c}$$\overline{b}c$$\overline{b}\overline{c}$$abc$$ab\overline{c}$$a\overline{b}c$$a\overline{b}\overline{c}$$\overline{a}bc$$\overline{a}b\overline{c}$$\overline{a}\overline{b}c$$\overline{a}\overline{b}\overline{c}$
Figure 43: Search graph for the Quine-McCluskey algorithm on 3 variables

Figure 44 shows the timing for our Quine-McCluskey run. Enumeration161616 We
still used the old enumeration scheme, i.e. we enumerated and checked also
relations in non-normal form.  of all relations on a $5$-element set and
checking all properties of each relation was done after $20$ seconds wall-
clock time. After that, levels $1$ to $5$ were completely checked within an
hour, but it took over a day to find the last law, on level $8$. The largest
share of run time was used in looping though all possible rectangles of a
level; level $n$ has $\left(\begin{array}[]{c}24\\\ n\end{array}\right)\cdot
2^{n}$ rectangles.

Note that the improved relation enumeration described in Sect. 6.1 would have
affected only the very first phase, which was completed in $20$ seconds,
anyway. However, it would have allowed for using a $6$-element universe set in
reasonable computation time, thereby avoiding the report of all non-laws in
Sect. 5, except 186 (Exm. 76), which needs $\geqslant 7$ elements in the
universe.

Somewhat unexpected, the $2^{25}$ relations inhabited no more than $495$ of
the $2^{24}$ possible combinations of properties. A listing of property
combinations by number of satisfying relations is provided in the ancillary
file propertyCombinationsByCount.txt.

Time |
---|---
(sec) |
0 | Counting property combinations
20 | Number of relations (Fig. 42)
21 | Prime implicants (Fig. 6–13)
| Level | On | Off | Don’t care
21 | 1 | 16776721 | 495 | 0
21 | 2 | 16776721 | 495 | 0
27 | 3 | 32063 | 495 | 16744658
85 | 4 | 161 | 495 | 16776560
575 | 5 | 32 | 495 | 16776689
3291 | 6 | 24 | 495 | 16776697
14591 | 7 | 4 | 495 | 16776717
49134 | 8 | 1 | 495 | 16776720
135925 | 9 | 0 | 495 | 16776721
Figure 44: Timing for Quine-McCluskey run

###  6.3 On finding “nice” laws

It is desirable to find a set of laws as “elegant” as possible. While
“elegance” is a matter of mathematicians’ taste and can hardly be rigorously
defined, some criteria for it are beyond doubt.

Each single law should be as general as possible. On a technical level, this
translates into two requirements:

  1. 1.

Each law should consist of as few literals as possible.

  2. 2.

Each law should use the sharpest predicates possible.

Considering sets of laws,

  1. 3.

some balance should be kept between conciseness and convenience.

We discuss criteria 1 to 3 in the following.

#### Criterion 1

Each law should consist of as few literals as possible.

That is, it should be obtained from a rectangle as large as possible. For
example, in a law $\lnot$ LfSerial $\lor$ $\lnot$ Asym $\lor$ Irrefl, the
first literal should be omitted.

This requirement is fulfilled, since our algorithm checks rectangles in order
of decreasing size, and for every reported rectangle prevents properly
contained rectangles from being reported also.

#### Criterion 2

Each law should use the sharpest predicates possible.

For example, both “LfUnique $\land$ RgUnique $\land$ SemiOrd2 $\Rightarrow$
ASym” (226) and “LfUnique $\land$ RgUnique $\land$ SemiOrd2 $\Rightarrow$
Irrefl” is a law, but the latter is redundant since it follows from the former
and the law “ASym $\Rightarrow$ Irrefl” (039).

In a naive approach to cope with this requirement, we ordered the properties
by extension set cardinality, see Fig. 42,171717 We used the figures from
column “Old” for that. Note that column “Pruned” would result in a slightly
different order. Also note that due to the nature of the pruning described in
Sect. 6.1, dual properties could have different extension set cardinalities.
and ensured that the algorithm checks rectangles in order of increasing
encoding. This way, 226 corresponds to the encoding 0x08000260, while its
weaker consequence would correspond to a larger encoding 0x08002060 and
therefore isn’t found by the algorithm.

However, when a predicate occurs negated in a law, this order doesn’t lead to
the desired result. For example, the law “LfUnique $\land$ Irrefl
$\Rightarrow$ AntiTrans” (127) is found, but its weaker consequence “LfUnique
$\land$ ASym $\Rightarrow$ AntiTrans” (109) was found before it; the former
and the latter corresponds to the encoding 0x002120 and 0x000320,
respectively.

The latter problem is caused by our too simple enumeration method. In
procedure qmc, we iterate in the loop on mask over all property sets of
cardinality given by the current level. For each such set we the iterate in
the loop on val over all assignments of negation symbols to the properties.
For example, referring to level 2 of Fig. 43, mask may take the values $ab$,
$ac$, $bc$, and for each value $xy$, val may take the values $xy$,
$x\overline{y}$, $\overline{x}y$, $\overline{x}\overline{y}$. However,
assuming that the order by extension set cardinalities is
$\overline{c}<\overline{b}<\overline{a}<a<b<c$, we should check the set $ac$
before the larger set $bc$ but after the smaller set $\overline{b}c$, which is
impossible with our simple enumeration method.

Figure 45: Partial rectangle order Induced by
$\overline{c}<\overline{b}<\overline{a}<a<b<c$

Given the cardinality of each property’s extension set (as in Fig. 42), a
partial order on the intersection set is induced; Fig. 45181818 Capital
letters denote negations.  gives an example for the facts about two-set
intersections inferrable191919 E.g. $a<b$ implies $ac<bc$ by monotonicity.
from the above order $\overline{c}<\overline{b}<\overline{a}<a<b<c$.

An improved approach should enumerate the rectangles in an order that is some
linearization of this inferred partial order. An efficient method to do this
is still to be found. One possibility might be to assign to each property a
weight, such that increasing extension set cardinalities correspond to
increasing weights, and to linearize the intersections in order of increasing
weight sums. Choosing powers of $2$ as weights will guarantee that all weight
sums are distinct; cf. the numbers in Fig. 45. However, sorting a list of all
$\left(\begin{array}[]{c}24\\\ n\end{array}\right)\cdot 2^{n}$ rectangles on
level $n$ would definitely not be efficient.

#### Criterion 3

Considering sets of laws, some balance should be kept between conciseness and
convenience.

For example, in a textbook about commutative groups, commutative variants of
the associativity axiom, like $(xy)z=(yz)x$, usually aren’t explicitly
mentioned, in order to keep the presentation concise. On the other hand, while
all theorems are redundant in the presence of an axiomatization, the book will
undoubtedly present some of them for convenience.

In our setting, we have a simple formal criterion about which laws to consider
redundant: those that follow from other laws solely by propositional logic.
For example, “CoRefl$\Rightarrow$LfEucl” (006), is not considered redundant
despite the triviality of its proof, since the latter needs to use Def. 1.3
and 1.16. In contrast,
“LfUnique$\land$IncTrans$\Rightarrow$ASym$\lor$LfSerial” (239) is redundant
since it follows from “LfUnique$\land$SemiOrd2$\Rightarrow$ASym$\lor$LfSerial”
(242) and “SemiOrd2$\Rightarrow$IncTrans” (071) by propositional inference
alone, without employing Def. 1.

rfiras039039039039039039– – 044044044044044044– – 046046046046046046– –
Figure 46: Redundant law suggestion example

As can be seen from the previous example, our algorithm doesn’t avoid
reporting redundant law suggestions. The reason for this is that it just
reports prime implicants in order of appearance, rather than selecting a
minimal covering subset of them. The latter technique is commonly employed in
proper Quine-McCluskey implementations, it is, however, NP-complete Feldman
[6, p.14]. The Karnaugh diagram in Fig. 46 illustrates the problem in a
simplified setting ($3$ properties only) along the example
“ASym$\Rightarrow$Irrefl” (039), “$\lnot$Refl$\lor\lnot$ASym” (044), and
“$\lnot$Refl$\lor\lnot$Irrefl” (046). Rectangle 044 is inspected before 046,
so the former is reported and its fields are set to don’t care. When the
latter is inspected, there is still one on field in it, so it is reported,
too. In contrast, a minimal covering subset approach wouldn’t report 044 since
all its fields are covered by the union of 039 and 046.

##  7 References

## References

  * Burghardt [2002] J. Burghardt, Axiomatization of Finite Algebras, in: Proc. KI 2002, no. 2479 in LNAI, Springer, ISBN 3-540-44185-9, 222–234, URL https://arxiv.org/abs/1403.7347, 2002\.
  * Quine [1952] W. Quine, The Problem of Simplifying Truth Functions, American Mathematical Monthly 59 (8) (1952) 521–531, URL https://www.jstor.org/stable/2308219.
  * McCluskey Jr. [1956] E. J. McCluskey Jr., Minimization of Boolean Functions, Bell System Technical Journal 35 (6) (1956) 1417–1444, URL https://archive.org/details/bstj35-6-1417.
  * Sen [1969] A. K. Sen, Quasi-Transitivity, Rational Choice and Collective Decisions, Review of Economic Studies 36 (3) (1969) 381–393, URL https://www.jstor.org/stable/2296434.
  * Sen [1971] A. K. Sen, Choice Functions and Revealed Preference, Review of Economic Studies 38 (3) (1971) 307–317, URL https://www.ihs.ac.at/publications/eco/visit_profs/blume/sen.pdf.
  * Feldman [2009] V. Feldman, Hardness of Approximate Two-Level Logic Minimization and PAC Learning with Membership Queries, Journal of Computer and System Sciences 75 (2009) 13–25, URL https://www.sciencedirect.com/science/article/pii/S002200000800069X/pdf?md5=0adeeefb189024849aaf393af8c98634&pid=1-s2.0-S002200000800069X-main.pdf.