axioms

Article
Local Influence for the Thin-Plate Spline Generalized
Linear Model
Germán Ibacache-Pulgar 1,2, * , Pablo Pacheco 3 , Orietta Nicolis 4 and Miguel Angel Uribe-Opazo 5

1 Institute of Statistics, Universidad de Valparaíso, Av. Gran Bretaña 1111, Valparaíso 2360102, Chile
2 Centro de Estudios Atmosféricos y Cambio Climático (CEACC), Universidad de Valparaíso,
Valparaíso 2360102, Chile
3 Dirección de Educación Virtual, Universidad de Playa Ancha, Avenida Guillermo González de Hontaneda
855, Playa Ancha, Valparaíso 2360072, Chile; ppachecog@upla.cl
4 Facultad de Ingenieria, Universidad Andres Bello, Calle Quillota 980, Viña del Mar 2520000, Chile;
orietta.nicolis@unab.cl
5 Centro de Ciências Exatas e Tecnológicas, Western Paraná State University (UNIOESTE),
Cascavel 85819-110, Paraná, Brazil; miguel.opazo@unioeste.br
* Correspondence: german.ibacache@uv.cl

Abstract: Thin-Plate Spline Generalized Linear Models (TPS-GLMs) are an extension of Semipara-
metric Generalized Linear Models (SGLMs), because they allow a smoothing spline to be extended to
two or more dimensions. This class of models allows modeling a set of data in which it is desired
to incorporate the non-linear joint effects of some covariates to explain the variability of a certain
variable of interest. In the spatial context, these models are quite useful, since they allow the effects
of locations to be included, both in trend and dispersion, using a smooth surface. In this work, we
extend the local influence technique for the TPS-GLM model in order to evaluate the sensitivity of
the maximum penalized likelihood estimators against small perturbations in the model and data.
We fit our model through a joint iterative process based on Fisher Scoring and weighted backfitting
algorithms. In addition, we obtained the normal curvature for the case-weight perturbation and
response variable additive perturbation schemes, in order to detect influential observations on the
model fit. Finally, two data sets from different areas (agronomy and environment) were used to
illustrate the methodology proposed here.

Citation: Ibacache-Pulgar, G.; Keywords: exponential family; smoothing spline; penalized likelihood function; weighted back-
Pacheco, P.; Nicolis, O.; Uribe-Opazo, fitting algorithm; diagnostics measures
M.A. Local Influence for the
Thin-Plate Spline Generalized Linear MSC: 62P12; 62J20; 62G05
Model. Axioms 2024, 13, 346. https://
doi.org/10.3390/axioms13060346

Academic Editor: Giovanni Nastasi

1. Introduction
Received: 19 April 2024 Thin-Plate Spline Generalized Linear Models (TPS-GLMs) represent an extension
Revised: 12 May 2024
of semiparametric generalized linear models (SGLMs) by enabling the application of
Accepted: 20 May 2024
smoothing splines in multiple dimensions. These models have the same characteristics
Published: 23 May 2024
of the generalized linear model (GLM), as described by McCullagh and Nelder [1]. Like
GLMs, TPS-GLMs can assume a variety of distribution families for the response variable.
They also allow for a non-linear relationship between the response variable’s mean and
Copyright: © 2024 by the authors.
the linear predictor via a link function, and they account for non constant variance in the
Licensee MDPI, Basel, Switzerland. data. Furthermore, the TPS-GLM allow modeling non-linear joint interaction effects due to
This article is an open access article some covariates, as well as the effects of coordinates in spatial data, making them a useful
distributed under the terms and tool to model dynamic pattern in different scientific areas, such as environment, agronomy,
conditions of the Creative Commons ecology, and so on. Some of the main works related to thin-plate spline technique are
Attribution (CC BY) license (https:// Duchon [2,3], Bookstein [4], and Chen et al. [5], while in the context of statistical modeling,
creativecommons.org/licenses/by/ Wahba [6], Green and Silverman [7], Wood [8], and Moraga et al. [9], can be mentioned,
4.0/). among others.

Axioms 2024, 13, 346. https://doi.org/10.3390/axioms13060346 https://www.mdpi.com/journal/axioms

Axioms 2024, 13, 346 2 of 20

However, it is well known that diagnostic analysis is a fundamental process in all

statistical modeling for any data set. This analysis allows us to validate the assumptions
established about the model in question and identify discrepant observations, and even-
tually influential ones on the fit of the model. One of the main diagnostic techniques
used in GLM and SGLM is local influence. In general, the idea of the local influence
technique introduced by Cook [10] is to evaluate the sensitivity of the MLEs when small
perturbations are introduced in the assumptions of the model or in the data, both in the
response variable and in the explanatory variables. This technique has the advantage,
regarding the case elimination technique, that it is not necessary to calculate the estimates
of the parameters for each case excluded. In our case, we are interested in developing
the local influence technique in the TPS-GLM, in order to detect observations that may
have a disproportionate influence on the estimators of both the parametric (regression
coefficient) and non-parametric (surface) part of the linear predictor. Such influence may
be due, for example, to the fact that each experimental unit contributes differently to the
model or that our variable of interest is exposed to a certain modification. In the context
of GLM and SGLM, there is empirical evidence that the maximum likelihood estimators
(MLEs) and maximum penalized likelihood estimators (PMLEs) are sensitive to this type of
situation, and therefore we believe that this sensitivity is also present in the estimators of
the TPS-GLM, in particular, in the surface estimator.
Various studies have expanded upon the technique of local influence within different
parametric models. Thomas and Cook [11] applied Cook’s method of local influence [10]
to generalized linear models to assess the impact of minor data perturbations. Ouwens
and Beger [12] obtained the normal curvature under a generalized linear model in order to
identify influential subjects and/or individual observations. Zhu and Lee [13] developed
the local influence technique for incomplete data, and extended such results to generalized
linear mixed models (see also Zhu and Lee [14] for further details). Espinheira et al. [15]
extended the local influence analysis to beta regression models considering various pertur-
bation scenarios. Rocha and Simas [16] and Ferrari et al. [17] derived the normal curvature
considering a beta regression model whose dispersion parameter varies according to the
effect of some covariates. Ferreira and Paula [18] developed the local influence approach to
partially linear Skew Normal models under different perturbation schemes, and Emami [19]
evaluated the sensitivity of Liu penalized least squares estimators using local influence
analysis. Most recently, Liu et al. [20] have reported the implementation of influence
diagnostics in AR time series models with Skew Normal (SK) distributions.
Within a semiparametric framework, Thomas [21] developed diagnostics for local
influence to assess the sensitivity of estimates for the smoothing parameter, which were
determined using the cross-validation criterion. Zhu and Lee [14] and Ibacache-Pulgar and
Paula [22] introduced measures of local influence to analyze the sensitivity of maximum
penalized likelihood estimates in normal and partially linear Student-t models, respectively.
Ibacache-Pulgar et al. [23,24] explored local influence curvature within elliptical semi-
parametric mixed models and symmetric semiparametric additive models. Subsequently,
ref. [25] and Ibacache-Pulgar and Reyes [26] further extended local influence measures to
normal and elliptical partially varying-coefficient models, respectively. Ibacache-Pulgar
et al. [27] developed the local influence method within the context of semiparametric addi-
tive beta regression models. Meanwhile, Cavieres et al. [28] calculated the normal curvature
to assess the sensitivity of estimators in a thin-plate spline model that incorporates skew
normal random errors. Jeldes et al. [29] applied the partially coefficient-varying model
with symmetric random errors to air pollution data from the cities of Santiago, Chile, and
Lima, Peru. In this context, they carried out an application of the local influence technique
to detect influential observations in the model fit. Saavedra-Nievas et al. [30] extended the
local influence technique for the spatio-temporal linear model under normal distribution
and with separable covariance. Recently, Sánchez et al. [31] obtained the normal curvature
for the varying-coefficient quantile regression model under log-symmetric distributions,
Axioms 2024, 13, 346 3 of 20

and presented an interesting application of such results to an environmental pollution

data set.
In this work, we extend the local influence approach in Thin-Plate Spline Generalized
Linear Model.
The contents are organized as follows: Section 2 introduces the thin-plate spline
generalized linear model. Section 3 details the method for obtaining maximum penalized
likelihood estimators and discusses some statistical inferential results. In Section 4, we
provide a detailed description of the local influence method and derives normal curvatures
for various perturbation schemes. In Section 5, the methodology is illustrated using two
datasets. The paper concludes with some final observations in Section 6.

2. The Thin-Plate Spline Generalized Linear Model (TPS-GLM)

In this section, we present the TPS-GLM and the penalized function to carry out the
process of estimating the parameters.

2.1. Statistical Model

Let {yi | i = 1, . . . , n} be a data set where each response variable yi follows a distribu-
tion from the exponential family with the following density function:
" #
yi θi − ψ ( θi )
f y (yi ; θi , ϕ) = exp + c ( yi , ϕ ) ,
ai ( ϕ )

where θi is the canonical form of the location parameter and depends on the mean µi .
The term ai (ϕ) represents a known function of the unknown dispersion parameter ϕ
(or a vector of unknown dispersion parameters). The function c depends on both the
dispersion parameter and the responses, while ψ is a known function, such that the mean
and variance of yi are given by: µi = E(yi ) = ∂ψ(θi )/∂θi and Var(yi ) = ai (ϕ) Vi , with
Vi = V(µi ) = ∂2 ψ(θi )/∂θi2 , respectively. The TPS-GLM is defined by Equation (1) and the
following systematic component:

g(µi ) = ηi = wi⊤ α + f (ti ), (1)

where wi is a ( p × 1) vector of covariables, α = (α1 , . . . , α p )⊤ corresponds to the vector of re-

gression coefficients, f (·) is unknown smooth arbitrary surface, and ti is a two-dimensional
covariates vector. To write the model given by Equation (1) in a matrix form, first con-
sider the one-to-one transformation of the vector f suggested by Green and Silverman [7],
stated as

f ( t1 )
 
.. T
f= .  = Eδ + T a ,
 

f (tn )

where a is a 3 × 1 vector with components ai , δ is a n × 1 vector with components δi , E is a

1
(n × n) matrix whose elements are given by Eij = 16π ∥ti − t j ∥2 log ∥ti − t j ∥2 , with Eii = 0
for each i, and T is a (3 × n) matrix defined as

1 1 ... 1
 
T= .
t1 t2 ... tn

Thus, the Model (1) can be written in a matrix form as

η = X β + Eδ,

T
T
where the regression matrix is structured as X = TT W , with W = w1⊤ , . . . , w⊤
n ,

T
and the vector of regression coefficients as β = aT α = β 1 , . . . , β p+3 , where


Axioms 2024, 13, 346 4 of 20

β j = a j (j = 1, 2, 3) and β j = α j−3 (j = 4, . . . , p + 3); see [9]. Note that this matrix

representation of the linear predictor allows us to treat the TPS-GLM as a semiparametric
generalized linear model, in which the term X β represents the parametric component and
Eδ the nonparametric component. One of the advantages of the TPS-GLM, apart from
being able to model both discrete and continuous variables that belong to the exponential
family, is its flexibility to model the non-linear joint effect of covariates through the surface
f present in the linear predictor η. In the context of spatial data, this models allows the
effect of coordinates to the incorporated into the modeling process. It is important to note
that when the surface f is not present in the linear predictor η, the model reduces to the
classical generalized linear model. However, if the vector t reduces to a scalar, t, the model
reduces to the semiparametric generalized linear model discussed, for instance, by Green
and Silverman [7].

2.2. Penalized Function

∗
Under the TPS-GLM, we have that θ = ( β⊤ , δ⊤ , ϕ)⊤ ⊆ R p , with p∗ = ( p + 3) + n + 1
parameters. Then, the log-likelihood function is given by
n
L(θ) = ∑ Li (θ), (2)
i =1

where
" #
yi θi − ψ ( θi )
Li (θ) = + c ( yi , ϕ ) .
ai ( ϕ )

To ensure the identifiability of the parameter vector α, we assume that f belongs to the
function space where all partial derivatives of total order m reside within the Hilbert space
L2 [ Ed ], the space of square-integrable functions on Euclidean d-space. Incorporating a
penalty function over f , we have that the penalized log-likelihood function can be expressed
as (see, for instance, Green and Silverman [7])
d
Lp (θ, λ f ) = L(θ) + λ∗f Jm (f) , (3)

where Jmd ( f ) is a penalty functional measuring the wiggliness of f , and λ∗ ( λ ) is a constant

f f
that depends on the smoothing parameter λ f ≥ 0. In general, a measure of the curvature
of f corresponds to its squared norm, ∥ f ∥, defined as
!2
d
∂m f
Z +∞ Z +∞
d m!
Jm (f) = ∥f∥ = ∑ υ1 ! . . . υ d ! −∞
...
−∞
α1 α
∂t1 . . . ∂td d
∏ dtj .
υ1 +...+υd =m j =1

For simplicity, in this work, we will consider the case in which d = 2, m = 2 and
g = g(t1 , t2 ). Consequently, the penalty function J22 ( f ) is expressed in the form
!2 !2 !2 )
∂2 f ∂2 f ∂2 f
Z Z
(
J22 ( f ) = +2 + dt1 dt2 ,
R2 ∂t21 ∂t1 ∂t2 ∂t22

and measures the rapid variation in f and the departure from local linearity. In this case, the
estimation of f leads to a natural thin-plate spline. According to Green and Silverman [7],
we may express the penalty functional as J22 ( f ) = δ T Eδ. Then, if we consider λ∗f = −λ f /2,
the penalized log-likelihood function (3) can be expressed as

λf T
Lp (θ, λ f ) = L(θ) − δ Eδ . (4)
2
Axioms 2024, 13, 346 5 of 20

The first term in the right-hand side of Equation (4) measures the goodness-of-fit,
while the second terms penalizes the roughness of f with a fixed parameter λ f . Selecting
appropriate parameters is crucial in the estimation process, as they determine the balance
between the goodness-of-fit and the smoothness (or regularity) of the estimated function. It
is important to emphasize that selecting appropriate parameters is crucial in the estimation
process because they control the trade-off between goodness-of-fit and the smoothness (or
regularity) of the estimated function. In this work, the smoothing parameter is selected
through the Akaike Criterion (AIC) based on the penalized log-likelihood function given in
Equation (3). More details of the method are given in Section 3.7.

3. Estimation and Inference

In this section, we discuss the problem of estimating the parameters under the TPS-
GLM. Specifically, we derive a weighted iterative process based on the backfitting algorithm
and estimate the variance–covariance matrix of our estimator from the penalized Fisher
information matrix (see Green [32] and Green and Silverman [7]). A brief discussion of the
smoothing parameter selection is also presented.

3.1. Penalized Score Function

First, we are going to assume that the function Lp (θ, λ f ) is regular in the sense that it
admits first and second partial derivatives with respect to the elements of the parameter
vector θ. To obtain the score function for β, we must calculate ∂Lpi (θ, λ f )/∂β j for i ∈
{1, . . . , n} and j ∈ {1, . . . , p + 2}. After performing some partial derivative operations, we
have that the score function for β can be written in matrix as follows:

β ∂Lp (θ, λ f )
Up (θ) = = X⊤ T
e (y − µ ) ,
∂β
h   i
e = diag (ai (ϕ))−1 ∂µi V −1 is a
where X is an (n × 3 + p) matrix whose ith row is xi⊤ , T ∂η i i
(n × 3 + p) matrix, with Vi = V (µi ) = ∂2 ψ(θi )/∂θi2 the variance function, ai (ϕ) a function
of ϕ, y = (y1 , ..., yn )⊤ and µ = (µ1 , ..., µn )⊤ are (n × 1) vectors.
Conversely, to derive the score function for δ, we need to compute ∂Lpi (θ, λ f )/∂δℓ for
i ∈ {1, . . . , n} and ℓ ∈ {1, . . . , n}. Again, after some algebraic operations, the score function
for δ can be written in matrix as follows:
∂Lp (θ, λ f )
Upδ (θ) = = E⊤ T
e (y − µ) − λ f Eδ ,
∂δ
where the matrix E is defined in Section 2.1. Finally, the score function for ϕ is given by

∂Lp (θ, λ f ) n n
= − ∑ (ai (ϕ))−2 {yi θi − ψ(θi )} + ∑ c′ (yi , ϕ),
ϕ
Up (θ) =
∂ϕ i =1 i =1

with c′ (yi , ϕ) = ∂c(yi , ϕ)/∂ϕ, for i ∈ {1, . . . , n}. Thus, the vector of penalized score func-
tions of θ can be expressed compactly as
β
 
Up (θ)
Up (θ) =  Upδ (θ)  .
 
ϕ
Up (θ)

Note that if the model under consideration only considers the parametric component
in the linear predictor, that is, the nonparametric component is omitted, the expressions of
the remaining score functions are reduced to those obtained under the classical generalized
linear model.
Axioms 2024, 13, 346 6 of 20

3.2. Penalized Hessian Matrix

To obtain the penalized Hessian matrix, we must compute the second-derivate of
Lp (θ, λ f ) with respect to each element of θ, that is, ∂2 Lp (θ, λ f )/∂θ j∗ θℓ∗ , for j∗ , ℓ∗ ∈ {1, . . . , p∗ }.
After some algebraic operations, we have that the diagonal elements (block matrices) of the
Hessian matrix are given by

ββ ∂2 Lp (θ, λ f )
Lp = = −X⊤ M∗ X ,
∂β∂α⊤
∂2 Lp (θ, λ f )
Lpδδ = ⊤
= −E⊤ M∗ E − λ f E and
∂δ∂δ
∂2 Lp (θ, λ f ) n n
∑ 2(ai (ϕ))−3 (yi θi − ψ(θi )) + ∑ c′′ (yi , ϕ)) ,
ϕϕ
Lp = =
∂ϕ2 i =1 i =1

where M∗ = diag1≤i≤n (ai (ϕ))−1 (∂µi /∂ηi )2 Vi−1 and c′′ (yi , ϕ) = ∂2 c(yi , ϕ)/∂ϕ2 , for 1 ≤
h i

i ≤ n. The elements outside the main diagonal of the Hessian matrix take the form

βδ ∂2 Lp (θ, λ f )
Lp = = −X⊤ M∗ E ,
∂β∂δ⊤
∂2 Lp (θ, λ f ) n
( )
βjϕ ∂µ
Lp = = − ∑ (ai (ϕ))−2 (yi − µi )Vi−1 i xij and
∂α j ∂ϕ i =1
∂ηi

∂2 Lp (θ, λ f ) n
( )
∂µ
= − ∑ (ai (ϕ)) −2
(yi − µi )Vi−1 i eiℓ
δ ϕ
Lpℓ = ,
∂δℓ ∂ϕ i =1
∂ηi

where xij denotes the (i, j)th element of the matrix X and eij denotes the (i, ℓ)th element of
the matrix E, for i ∈ {1, . . . , n}, j ∈ {1, . . . , p + 2} and ℓ ∈ {1, . . . , n}. Thus, the penalized
Hessian matrix can be represented as
ββ βδ βϕ
 
Lp Lp Lp
 βδ⊤ δϕ
Lp (θ) =  .

 Lp Lpδδ Lp 
βϕ⊤ δϕ⊤ ϕϕ
Lp Lp Lp

It is noteworthy that this matrix simplifies to the Hessian matrix used in generalized
linear models when the nonparametric component is absent. The primary application of
this matrix lies in the normal curvature, which is essential for developing the local influence
technique. This will be discussed in the following section.

3.3. Penalized Expected Information Matrix

By taking the expectation of the matrix −Lp (θ), we derive the penalized expected
information matrix, which is of dimension ( p∗ × p∗ ), as follows:

∂2 Lp (θ, λ)
" #
Jp (θ) = −E .
∂θ∂θT

This matrix assumes the following diagonal structure in blocks:

!
βδ
Jp (θ) 0
Jp (θ) = ϕϕ ,
0 Jp (θ)
Axioms 2024, 13, 346 7 of 20

where

X⊤ M∗ X X⊤ M∗ E
 
βδ
Jp (θ) =
E M X E M∗ E + λ f E
⊤ ∗ ⊤

and
n n
∑ −2(ai (ϕ))−3 (µi θi − ψ(θi )) − ∑ E(c′′ (yi , ϕ)),
ϕϕ
Jp (θ) =
i =1 i =1

with c′′ (yi , ϕ) = ∂2 c(yi , ϕ)/∂ϕ2 for i ∈ {1, . . . , n}.

3.4. Derivation of the Iterative Process

The value of θ that maximizes Lp (θ, λ f ), called maximum penalized likelihood es-
timate (MPLE) and denoted by b θ, is carried out by solving the corresponding estima-

⊤  ⊤
tion equations. Let θ = θ1 θ2 , where θ1 = β⊤ δ⊤
⊤ and θ2 = ϕ. In addition,
 ⊤
⊤ ⊤
consider the partition of the score function vector Up (θ) = U1p (θ) U2p (θ) , where
 ⊤ 
⊤ ϕ
U1p (θ) = Upβ (θ) Upδ (θ) and U2p (θ) = Up (θ). In order to estimate θ based on penali-
⊤

zed likelihood function given by Equation (4), we have to solve the equations

U1p (θ) = 0


U2p (θ) = 0.

These estimating equations are nonlinear, and necessitate an iterative approach for
their solution. An alternative frequently proposed in the context of generalized linear
models is the Fisher scoring algorithm (Nelder and Wedderburn, [33]), considering the fact
that in some situations the matrix −Lp (θ) can be non-positive definite. Then, the algorithm
for estimating θ1 , with ϕ fixed, is given by
old
θnew = θold −1 old 1
βδ
1 1 + (Jp ( θ) ) Up (θ),

which is equivalent to solving the matrix equation

Sold βnew
!   old old 
I β E S0 z
new = , (5)
Sold
δ X I δ Sold
1 z
old

where zold = (y − µold ) + ηold , with Sold

ϑ defined as

old old

 (X⊤ M∗ X ) −1 X⊤ M∗
 ϑ=β
Sϑ =
 (E⊤ M∗old E + λ E)−1 E⊤ M∗old

ϑ=δ.
f

Consequently, the weighted back-fitting (Gauss–Seidel) iterations for simultaneously

updating β and δ are given by

βnew = Sold zold − Eδold , (6)

β
new
= Sold zold − Xβold , (7)


δ δ

It is crucial to note that the system of Equations (5) is consistent, and the back-fitting
algorithm converges to a solution for any initial values, provided that the weight matrix
Axioms 2024, 13, 346 8 of 20

M∗ is symmetric and positive definite. Additionally, if the parametric component wi⊤ β is

absent in the linear predictor, the estimator of δ is given by:

δnew = Sold
δ z
old
.

The MPLE of the dispersion parameter, θb2 = ϕ

b, can be determined through the
following iterative procedure:
old
θ2new = θ2old − (Jp (θ)−1 )old U2p (θ) .
ϕϕ

Summarizing, each iteration of the Fisher scoring algorithm updates β and δ using
Equations (6) and (7), and evaluating matrices Sϑ and M∗ at the MPLE of θ obtained in the
previous iteration, that is, θold , until convergence is obtained. The joint iterative process
that resolves Up (θ) = 0 is presented below.

3.5. Estimation of Surface

To obtain the MPLE of f, we must consider its one-to-one representation given in
Equation (2) and MPLE obtained from the iterative process described above. Indeed, we
have that bf can be obtained as

δ + TT b
bf = Eb a, ,

where b δ and ba are the MPLE of δ and b

a, respectively. Note that vector ba corresponds to the
first three elements of vector β.
b Consequently, bf is a natural thin-plate spline. Details of the
conditions that guarantee this result are given, for example, in Green and Silverman [7]
and Wood [34].

3.6. Approximate Standard Errors

In this study, we propose approximating the variance–covariance matrix of bθ by using
the inverse of the penalized Fisher information matrix. Specifically, we have that

Cov
d (bθ) ≈ Jp−1 (θ) bθ .

If we are interested in drawing inferences for β, the approximate variance–covariance

matrix can be estimated by using the corresponding block-diagonal matrix obtained from
Jp−1 (θ), similarly for f and ϕ.

3.7. On Degrees of Freedom and Smoothing Parameter

For the TPS-GLM, the degree of freedom (d f ) associated with the smooth surface is
given by (see, for instance Green and Silverman [7])

d f (λ f ) = tr(E⊤ Sδ ),

which approximately represents the number of effective parameters used in the modeling
process to estimate the smooth surface f .
Regarding the selection of the smoothing parameter, we propose to use the Akaike
Information Criterion (AIC) (see, for instance, [24,35]), defined as

AIC (λ f ) = −2L p (θ, λ f ) + 2[1 + p + d f (λ f )] ,

θ
b

where L p (θ, λ f ) denote the penalized likelihood function evaluated at MPLE of θ, and p
denote the number of parameters in β. As usual, the idea is to select the value of λ f that
minimizes AIC (λ f ).
Axioms 2024, 13, 346 9 of 20

4. Local Influence
In this section, we extend the local influence technique to evaluate the sensitivity of
the MPLE under the TPS-GLM. Specifically, we present some theoretical aspects of the
method and, subsequently, we derive the normal curvature for three perturbation schemes.

4.1. Local Influence Analysis

Consider ω = (ω1 , . . . , ωn )⊤ , an n × 1 vector of perturbations restricted to some
open subset Ω ⊂ Rn . Let L p (θ, λ f | ω) denote the logarithm of the perturbed penalized
likelihood function. Assume there exists a vector of non-perturbation ω0 ∈ Ω, such that
L p (θ, λ f | ω0 ) = L p (θ, λ f ). To evaluate the influence of small perturbations on the MPL
estimate θ, b we can consider the penalized likelihood displacement given by:

2 [ L p (b θω , λ f )] ≥ 0,
θ, λ f ) − L p (b

where θ bω is the MPL estimate under L p (θ, λ f | ω). The measure LD(ω) is useful for
assessing the distance between θ
b and θ
bω . Cook [10] suggested examining the local behavior
of LD(ω) around ω0 . The procedure involves selecting a unit direction d ∈ Ω, with ∥d∥ = 1,
and then plotting LD(ω0 + ad) against a, where a ∈ R. This plot, called the lifted line, can be
characterized by considering the normal curvature Cd (θ) around a = 0. The suggestion is to
assume the direction d = dmax corresponding to the largest curvature Cdmax (θ). The index
plot of dmax can identify those cases that, under small perturbations, have a significant
potential influence on LD(ω). According to Cook [10], the normal curvature in the unit
direction d is expressed as

Cd (θ) = −2{d⊤ ∆p⊤ Lp−1 ∆p d},

with

∂2 L p (θ, λ f ) ∂2 L p (θ, λ f ω)
Lp = and ∆p = .
∂θ∂θ⊤ θ=θ
b ∂θ∂ω⊤ θ=θ,
b ω=ω0

Note that −Lp represents the penalized observed information matrix evaluated at
b (see Section 3.2), and ∆p is the penalized perturbation matrix evaluated at θ
θ b and ω0 .
It is essential to highlight that Cd (θ) denotes the local influence on the estimate θ b after
perturbing the model or data. Escobar and Meeker [36] suggested examining the normal
curvature in the direction d = ei , where ei is an n × 1 vector with a one at the ith position
and zeros elsewhere. Consequently, the normal curvature, referred to as the total local
influence of the ith case, takes the form Cei (θ) = 2|cii | for i ∈ {1, . . . , n}, where cii is the ith
principal diagonal element of the matrix C = ∆p⊤ Lp−1 ∆p .

4.2. Derivation of the Normal Curvature

Typically, the perturbation schemes used in the analysis of local influence are de-
termined by the structure of the model under consideration, as discussed by Billor and
Loynes [37]. These schemes can generally be divided into two main categories: perturba-
tions to the model (to examine changes in assumptions) or perturbations to the data. For
instance, we might consider perturbing the response variable or the explanatory variables.
The motivation for employing these perturbation schemes often includes addressing issues
such as the presence of outliers or the occurrence of measurement errors in the data. Subse-
quently, we will present the formulas for the matrix ∆p for various perturbation schemes.
Consider the weights assigned to the observations in the penalized log-likelihood
function, given by:
n λf ⊤
L p (θ, λ f ω) = L(θ ω) − ∑ δ Eδ ,
i =1
2
Axioms 2024, 13, 346 10 of 20

where L(θ ω ) = ∑in=1 ωi Li (θ), ω = (ω1 , . . . , ωn )⊤ is the vector of weights, with 0 ≤ ωi ≤

1. In this case, the vector of no perturbation is given by ω0 = 1(n×1) . Differentiating
L p (θ, λ f ω) with respect to the elements of θ and ω T , we have that the matrix ∆p takes
the form  ⊤
X Dτ


∆ p =  E⊤ Dτ  ,
b⊤
u

where the matrix Dτ = diag1≤i≤n τi and u = (u1 , . . . , un )⊤ , with τi = (ai (ϕ))−1 (yi −


′−1 ′−1 ′
∂ψ(h(ηi ))/∂h(ηi ))∂h(ηi )/∂ηi , h(ηi ) = ψ (ηi ), ψ (·) denotes the inverse function of ψ (·),
ui = −(ai (ϕ))−2 (yi h(ηi ) − ψ(h(ηi )) + c′ (yi , ϕ)ein
⊤ , and e a vector with 1 at the ith position
in
and zero elsewhere.
To perturb the response variable values, we consider yiω = yi + ωi for i ∈ {1, . . . , n},
where ω = (ω1 , . . . , ωn )⊤ is the vector of perturbations. The vector of no perturbation
is ω = (0, . . . , 0)⊤ . The perturbed penalized log-likelihood function is constructed from
Equation (3) with yi replaced by yiω , as follows:
n λf ⊤
L p (θ, λ g ω) = L(θ ω) − ∑ δ Eδ ,
i =1
2

where L(·) is defined in Equation (2) with yiω replacing yi . By differentiating L p (θ, λ | ω)
with respect to the elements of θ and ωi , and after some algebraic manipulation, we obtain:

X⊤ Dc
 

∆ p =  E⊤ Dc  ,
b⊤
d

where the matrix Dc = diag1≤i≤n ci and d = (d1 , . . . , dn )⊤ , with ci = ∂h(ηi )/∂ηi and

di = −(ai (ϕ))−2 (h(ηi )ein

⊤ + c′ ( y , ϕ ) /∂ω ), with e denoting a vector with 1 at the ith
iω i in
position and zero elsewhere.

5. Applications
In this section, we show the applicability of the TPS-GLM and the local influence
method with two real data applications. The model estimation and diagnostics have been
implemented using MatLab 9.13.0 (R2022b) software [38] (the developed code is available
on request by the authors).

5.1. Wypych Data

The first dataset we use to illustrate the applicability of the TPS-GLM consists of
83 sample points within a 46.6-hectare agricultural area in Wypych, located at latitude
24°50′ 24′′ S and longitude 53°36′ 36′′ W, with an average altitude ranging from 589 to 660 m.
The data were collected during the 2006/2007 agricultural year in the western region
of Paraná State, Brazil (see [39], Appendix 4). The soil is classified as Dystroferric Red
Latosol with a clayey texture. The region’s climate is mesothermal, super-humid temperate,
classified as Cfa according to (Köeppen), with a mean annual temperature of 21 °C. The
83 georeferenced points were determined by a regular grid of 75 × 75 m using a global
positioning system (GPS). The collected variables were as follows:
• Soya: average of soybean yield (t/ha).
• Height: average height (cm)of plants at the end of the production process.
• Pods: average number of pods.
• Lat: latitude (UTM).
• Long: longitude (UTM).
The original objective was to investigate the spatial variability of soybean yield (Soya)
in the studied area based on the covariates: average plant height, average number of pods
Axioms 2024, 13, 346 11 of 20

per plant, latitude, and longitude. Figure 1 shows the scatterplots between the response
variable Soya and the explanatory variables Height and Pods. In addition, the plot of the
response variable against the coordinates is shown. Clearly, from Figure 1a,b, it can be seen
that the explanatory variables Height (X2) and Pods (X3) are linearly related to the response
variable Soya (Z). The spatial effect given by the coordinates (X,Y) will be incorporated into
the model through a smooth surface.

(a) (b)
6 6

5.5 5.5

5 5

4.5 4.5

4 4
Z

Z
3.5 3.5

3 3

2.5 2.5

2 2

1.5 1.5
30 35 40 45 50 55 60 65 20 25 30 35 40 45 50 55 60 65
X2 X3

(c)

4
Z

1
2.37
2.368 7.2512
7.251
×10 5 2.366
7.2508 ×10 6
2.364
7.2506
X 2.362 7.2504 Y

Figure 1. Scatter plots: Soya versus Height (a), Soya versus Pods (b), and Soya versus coordinates in
UTM (c).

5.1.1. Fitting the TPS-GLM

Based on the above analysis, we propose the TPS-GLM, introduced in Section 2, to
model the trends present in the Wypych data. Specifically, we are going to assume that the
response variable Soya belongs to the Gaussian family, and that the link function is the
identity. Therefore, the model is expressed as follows:

g(µi ) = µi = β 0 + β 1 Heighti + β 2 Podsi + f (Lati , Longi ) i ∈ {1, . . . , 83},

where β = ( β 0 , β 1 , β 2 )⊤ correspond to the regression coefficients associated with the

parametric component of the model, and f (·) is a smooth surface. Table 1 lists the MPLE of
β. The respective asymptotic standard errors are presented in parentheses.
Axioms 2024, 13, 346 12 of 20

Table 1. MPLEs with their standard errors (within parenthesis), AIC and R2 (Adj).

Model
Parameters Gaussian Linear TPS-GLM
β0 1.1921 (0.672) 0.497 (0.751)
β1 0.0116 (0.0128) 0.032 (0.015)
β2 0.0339 (0.0079) 0.030 (0.008)
AIC 149.99 139.9992
R2 (Adj) 0.168 0.315

The value of the smoothing parameter λ f was selected in such a way that the AIC
criterion was minimized. The adjusted determinant coefficients (R2 (Adj)) are evaluated
for assessing the goodness-of-fit of the two models. It is important to note that our model
have a lower AIC and an higher R2 (Adj), compared to the multiple regression model that
does not consider the spatial effect. Figure 2a shows the QQ-plot for the standardized
residuals, whose adjustment to the Gaussian TPS-GLM seems to be reasonable. However,
the presence of some atypical observations is observed in one of the tails of the distribution.
Figure 2b displays the scatter plot between the observed values, Soya, and their estimated
values, Soya.
[ Considering the trend of the points, we conclude that the estimates are good,
since they generate consistent adjusted values of the response variable.

(a)
(b)
1.5
2.0
1.5

1
Empirical Quantiles

1.0

0.5
0.5
0.0

0
−0.5

-0.5
−1.0

−2 −1 0 1 2
-1
1.5 2 2.5 3 3.5 4 4.5 5 5.5 6
Norm Quantiles

Figure 2. QQ-plots of the standardized residuals for the TPS-GLM with its confidence interval
(dashed lines) (a) and scatterplot between Soya and Soya
[ (b), under model fitted to Wypych data.

5.1.2. Diagnostic Analysis

To identify potentially influential observations on the MPLE under the fitted Gaussian
TPS-GLM for the Wypych data, we present several index plots of Bi = Bei (γ) for γ = β, δ.
Figure 3 shows the index plot Bi for the case-weight perturbation scheme under the fitted
model. Figure 3a reveal that the observations #6, #61, #69 and #71 are more influential on
b whereas the observations #6, #66, #61 and #38 are more influential on b
β, δ; see Figure 3b.
When we perturb the response variable additively, we have that the observations #80, #32,
#75 and #88 are more influential on β;b see Figure 4a. Regarding b
δ, observations #3, #42 and
#80 appear as slightly influential as seem in Figure 4b.
Axioms 2024, 13, 346 13 of 20

(a) (b)
1 0.6

0.8 6
6

0.4
0.6

Bi

Bi
66
61
0.4 38
61 0.2

69 71
0.2

0 0
0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80
Index Index

Figure 3. Index plots of Bi for assessing local influence on β

b (a) and b
δ (b), considering case-weight
perturbation.

(a) (b)
0.6 0.4

0.4
80
80
Bi

Bi
32 0.2
3 42
88 75

0.2

0 0
0 10 20 30 40 50 60 70 80 0 10 20 30 40 50 60 70 80
Index Index

Figure 4. Index plots of Bi for assessing local influence on β

b (a) and b
δ (b), considering response
variable additive perturbation.

We conclude that the maximum penalized likelihood estimates (MPLE) of the regres-
sion coefficients and the smooth surface exhibit sensitivity to modifications made to the
data or the model. This analysis has shown that observations identified as influential for the
parametric component do not necessarily exert influence on the non-parametric component,
and vice versa. For instance, under the case-weight perturbation scheme, observations
#69 and #71 were detected as influential for the parametric component, but not for the
nonparametric component.

5.1.3. Confirmatory Analysis

Table 2 displays the relative changes experienced by each element in the vector of
regression coefficients. In this analysis, we only consider the three most influential observa-
tions under the case-weight perturbation scheme. As can be seen in this table, observations
#6, #61 and #69 generate significant changes in the estimates. Still, no relevant inferential
changes were noted. However, the AIC and &R2 (Adj) present some differences once the
above observations are dropped.
Axioms 2024, 13, 346 14 of 20

Table 2. Relative changes (RCs) (in %) in the MPL estimates of β j in cases-weight perturbation
under the TPS-GLM. The last two columns indicate the AIC and R2 (Adj) of the model with dropped
observations.

Parameters and Relatives Changes

Dropped Obs. β0 β1 β2 RCβ0 RCβ1 RCβ2 AIC R2 (Adj)
6 1.696 0.009 0.023 122.59 63.93 27.52 125.50 0.267
61 0.987 0.022 0.027 29.47 15.19 12.23 131.96 0.356
69 0.432 0.035 0.028 43.33 35.04 10.85 138.07 0.326
6-61 1.996 0.002 0.022 161.89 92.06 27.65 116.20 0.218
6-69 1.481 0.014 0.023 94.35 45.61 26.09 124.52 0.268
61-69 0.704 0.029 0.028 7.59 9.58 10.93 131.03 0.358
6-61-69 1.868 0.005 0.023 145.16 81.03 26.90 115.83 0.310

On the other hand, Table 3 shows the relative changes in the vector of regression
coefficients under the additive perturbation scheme of the response variable. Here, we
consider the four most influential observations. As can be seen in the table, observations
#32, #69, #75 and #80 generate important relative changes in the estimates of the parametric
component of the model. However, no significant inferential changes were observed. About
the AIC and &R2 (Adj), there are not evident differences.

Table 3. Relative changes (RCs) (in %) in the MPL estimates of β j in response variable perturbation
under the TPS-GLM. The last two columns indicate the AIC and R2 (Adj) of the model with dropped
observations.

Parameters and Relatives Changes

Dropped Obs. β0 β1 β2 RCβ0 RCβ1 RCβ2 AIC R2 (Adj)
32 0.701 0.027 0.029 8.000 4.62 5.23 138.73 0.319
69 0.432 0.034 0.027 43.33 29.0 12.22 138.07 0.326
75 0.760 0.028 0.027 0.24 6.22 12.65 139.44 0.307
80 0.699 0.028 0.029 8.21 7.86 8.16 139.01 0.311
32-69 0.382 0.035 0.030 49.87 32.82 4.00 136.81 0.33
32-75 0.700 0.027 0.030 8.17 4.12 4.90 138.11 0.319
32-80 0.621 0.028 0.031 18.49 6.34 0.16 137.63 0.316
69-75 0.430 0.035 0.028 43.61 34.96 10.96 137.53 0.318
69-80 0.333 0.036 0.029 56.33 38.32 5.39 136.99 0.322
75-80 0.695 0.028 0.029 8.85 7.25 7.90 138.39 0.302
32-69-75 0.381 0.035 0.030 49.98 32.33 3.74 136.22 0.322
32-75-80 0.621 0.028 0.031 18.54 5.23 0.96 136.93 0.308
69-75-80 0.333 0.036 0.029 56.26 37.40 5.09 136.42 0.314
32-69-75-80 0.271 0.035 0.032 64.47 30.15 3.22 134.96 0.320

5.2. Ozone Concentration Data

For our analysis, we utilize data from a study examining the relationship between
atmospheric ozone concentration (O3) and various meteorological variables in the Los
Angeles Basin for a sample of 330 days in 1976. The data were initially presented by
Breiman and Friedman [40], and are available for download from various public reposito-
ries. Although the dataset includes several variables, in this application, we will consider
only three explanatory variables, which are detailed in the following.
• O3: daily maximum one-hour average ozone concentration in Upland, CA, measured
in parts per million (ppm).
• Temp: Sandburg Air Base temperature, in Celsius.
• Vis: visibility, in miles.
• Day: calendar day.
Figure 5 contains the dispersion graphs between the outcome variable (log(O3)) and
each one of the explanatory variables Temp, Vis and Day.
Axioms 0 346
2024,1,13,
Axioms2024, 15
15of
of20
20

(a) (b)

4 4

3 3

log(O3)

log(O3)
2 2

1 1

0 0
100 400
80 400 300 400
60 300 200 300
200 200
40 100
100 100
Temp 20 0 Day Vis 0 0 Day

Figure5.5.3D
Figure 3Dplots
plotsbetween
betweenthe
theresponse
responsevariable
variableand
andthe
theexplanatory
explanatoryvariables:
variables:logarithm
logarithmof ofozone
ozone
data
data versus temperature and day variables (a), and logarithm of ozone data versus visibility andday
versus temperature and day variables (a), and logarithm of ozone data versus visibility and day
variables
variables(b).
(b).

Figure
Figure5a 5ashows
showsaacurved
curvedsurface
surfacein inthe
therelationship
relationshipbetween
betweenthethevariable
variablelog(O3)
log(O3)
and
and the joint effect of the explanatory variables Temp and Day, whereas therelationship
the joint effect of the explanatory variables Temp and Day, whereas the relationship
between
betweenlog(O3)
log(O3)and andthethejoint
jointeffect
effectof ofthe
theexplanatory
explanatoryvariable
variableVis
Visand
andDay
Dayshows
showsless
less
curve;
curve; see Figure 5b. This graphical analysis recommends the inclusion in the modelof
see Figure 5b. This graphical analysis recommends the inclusion in the model ofaa
nonparametric
nonparametriccomponent,
component,specifically
specificallyaasurface,
surface,that
thatcan
canexplain
explainthe
therelationship
relationshipbetween
between
log(O3)
log(O3)andandthe
thecombined
combinedeffecteffectof
ofthe
theexplanatory
explanatoryvariables
variablesTemp
TempandandDay.
Day.For
Forsimplicity,
simplicity,
in
inthis
thiswork,
work,we wewill
willinclude
includethetheeffect
effectofofthe
theexplanatory
explanatoryvariable
variableVis
Visininaalinear
linearform.
form.To To
begin
beginour
ouranalysis,
analysis,we wearearegoing
goingtotoconsider
considerthethefit
fitof
ofaaGLM
GLMassuming
assumingthatthatthe
thevariable
variableofof
interest O3 is Poisson distributed with mean and logarithmic link function.
interest O3 is Poisson distributed with mean µi i and logarithmic link function. Different
µ Different
structures
structuresof ofthe
thelinear
linearpredictor
predictorfor forthe
theexplanatory
explanatoryvariables
variablesVis,
Vis,Temp,
Temp,andandDay
Daywill
willbe
be
considered (see Table
considered (see Table 4). 4).

Table4.4. Four
Table Fourstructures
structuresof
ofthe
thelinear
linearpredictor
predictorfor
forthe
theexplanatory
explanatoryvariables
variablesVis,
Vis,Temp,
Temp,and
andDay,
Day,
assuming
assumingthatthatthe
theresponse
responsevariable
variablelog(O3)
log(O3)follows
followsaaPPOISSON
OISSON((µµi )i )distribution.
distribution.

Model
Model gg((µµi i))==log
log((µµi i))
II ββ00+ +ββ11Vis +ββ22Temp
Visi i+ +ββ33Day
Tempi i+ Dayi i
IIII ββ00+ + β11 Visi i + β22 Tempi i + ff((Day
β Vis + β Temp + Dayi i))
III
III ββ00+ +ββ11Vis +ββ22Temp
Visi i+ +ββ33Day
Tempi i+ Dayi i+ +ββ44Temp
Tempi i x×Day
Dayi i
IV
IV ββ00+ Vis + f ( Temp , Day )
+ β11 Visi i + f (Tempi i , Dayi i )
β

For
For Model
Model I,I, we we consider
consider onlyonly the
the individual
individual effects
effects of
of the
the explanatory
explanatory variables
variables
Vis,
Vis, Temp and Day. Note that all these effects were incorporated in a linear form
Temp and Day. Note that all these effects were incorporated in a linear form in in
the
the systematic component of the model. For Model II, we consider the inclusion of
systematic component of the model. For Model II, we consider the inclusion of aa
nonparametric
nonparametricterm termto tomodel
modelthe thenonlinear
nonlineareffects
effectsof
ofthe
theexplanatory
explanatoryvariable
variableDay;
Day;see see
Ibacache et al. [41]. Model III considers a systematic component that contains
Ibacache et al. [41]. Model III considers a systematic component that contains the individual the individual
effects
effectsof ofthetheexplanatory
explanatoryvariables
variablesVis,
Vis,Temp
TempandandDay,
Day,ininaddition
additionto tothe
theincorporation
incorporationof of
the
theinteraction
interactioneffect effectbetween
betweenthe theexplanatory
explanatoryvariables
variablesTemp
Tempand andDay.
Day.Here,
Here,thetheinteraction
interaction
effect
effectisisintroduced
introducedlinearly
linearlyin inthe
themodel.
model.Model
ModelIV IVcorresponds
correspondsto toaaTPS-GLM
TPS-GLMwhere wherethe the
joint
jointeffect
effectof ofthe
theTemp
Tempand andDayDayexplanatory
explanatoryvariables
variablesisisincluded
includednonlinearly
nonlinearlyby byusing
using
smooth
smoothsurface.
surface.Table
Table55contains
containsthetheML
MLandandMPLMPLestimates
estimatesassociated
associatedwithwiththetheparametric
parametric
component for the four fitted
component for the four fitted models. models.
ItItisisimportant
importantto tonote
notethat
thatboth
boththe
theindividual
individualandandinteraction
interactioneffects
effectsare
arestatistically
statistically
significant,
significant,as asthe
thecorresponding
correspondingp-values p-values(not
(notshown
shownhere)
here)are
areless
lessthan
than0.05.
0.05.Additionally,
Additionally,
the
theestimates
estimatesof of ββ00 are
aresimilar
similaracross
acrossthe
thefour
fourmodels,
models,whereas
whereasthe theestimates
estimatesof of ββ11 vary
vary
considerably,
considerably, particularly
particularly in in Model
Model IV.IV. Concerning
Concerning the the associated
associated standard
standard errors,
errors, all
all
Axioms 2024, 13, 346 16 of 20

the estimators exhibit small values. The last two rows of Table 5 display the Akaike
Information Criterion (AIC) and R2 values, respectively. It is evident that the TPS-GLM,
with AIC(λ f ) = 1777.705, provides the best fit to the Ozone data, followed by Model II
with an AIC of 1806.837. This is corroborated by the QQ-plots in Figure 6, specifically
Figure 6b,d. Furthermore, the R2 value associated with our model is higher than those
of Models I, II, and III. The smoothing parameter λ f was chosen such that the effective
degrees of freedom were approximately 7. Figure 7 illustrates the 3D plot of the adjusted
log(O3) against the explanatory variables Temp and Day, showing an adequate fit of the
TPS-GLM.

Table 5. AIC, R2 (Adj), ML and MPL estimates for all four fitted models to the Ozone data.

Parameters I II III IV
β0 0.577 (0.104) 0.478 (0.142) 0.787 (0.198) 2.507 (0.040)
β1 −0.002 (0.0003) −0.002 (0.0003) −0.002 (0.0003) −0.002 (0.0003)
β2 0.035 (0.001) 0.033 (0.002) 0.032 (0.003) -
β3 −0.001 (0.002) - −0.002 (0.001) -
β4 - - 0.00002 (0.00002) -
AIC 1887.312 1806.837 1887.757 1789.92
R2 (Adj) 0.673 0.715 0.670 0.728

(a) (b)

4
4

3
Empirical Quantiles

Empirical Quantiles

2
2

1
0

0
−1
−2

−2
−3

−3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3

Norm Quantiles Norm Quantiles

(c) (d)
3
4

2
Empirical Quantiles

Empirical Quantiles
2

1
0
0

−1
−2

−2

−3 −2 −1 0 1 2 3 −3 −2 −1 0 1 2 3

Norm Quantiles Norm Quantiles

Figure 6. QQ-plot of the standardized residuals for the models described in Table 5: Model I (a),
Model II (b), Model III (c) and Model IV (d).
Axioms 2024, 13, 346 17 of 20

(b)

log(µ ) 1

–1
100
80 400
60 300
200
40
100
Temp 20 0 Day

Figure 7. 3D plot between log

\ (µ) and explanatory variables Temp and Day.

5.2.1. Diagnostic Analysis

To identify potentially influential observations on the MPL estimators under the fitted
TPS-GLM for the Ozone data, we present some index plots of Bi = Bei (γ), for γ = β, δ.
Figure 8 shows the index plot of Bi for the case-weight perturbation scheme under the fitted
model. In Figure 8a,b, note that observations #167, #220, #168, and #177 are more influential
on β
b and bδ, respectively. By perturbing the response variable additively, it becomes clear
that observations #125, #175, #218, #219, and #221 are more influential on β b and b δ; see
Figure 9a and 9b, respectively.

(a) (b)
0.6 0.8

0.6 220
167
0.4 167 220
Bi

Bi

0.4
168 177
168177
0.2

0.2

0 0
0 50 100 150 200 250 300 0 50 100 150 200 250 300
Index Index

Figure 8. Index plots of Bi for assessing local influence on β

b (a) and b
δ (b), considering case-weight
perturbation under model fitted to Ozone data.

From the local influence analysis, we conclude that the MPLE of the regression co-
efficients and the smooth surface are sensitive to perturbations in the data or the model.
Furthermore, this analysis revealed that observations identified as influential for the para-
metric component are also influential for the nonparametric component, and vice versa. For
example, under the case-weight perturbation scheme, observations #167, #220, #168, and
#177 were found to be influential for both the parametric and nonparametric components.
Axioms 2024, 13, 346 18 of 20

(a) (b)
0.4 0.4

219
175
220221
125 167
218
Bi

Bi
0.2 219 0.2
125 175
218221

0 0
0 50 100 150 200 250 300 0 50 100 150 200 250 300
Index Index

Figure 9. Index plots of Bi for assessing local influence on β

b (a) and b
δ (b), considering response
variable additive perturbation.

5.2.2. Confirmatory Analysis

To investigate the impact on model inference when influential potentially observations
detected in the diagnostic analysis are removed, we present the relative changes (RCs)
in the MPL estimate of β j for j ∈ {1, 2} after removing the influential observations from
ξb−ξb( I )
the dataset (%). The RC is defined as RCξ = × 100%, where ξb( I ) denotes the MPL
ξb
estimate of ξ, with ξ = β j , after the corresponding observation(s) are removed according
to set I. Table 6 presents the RCs in the regression coefficient estimates after removing the
observations identified as potentially influential for the parametric component of the model.

Table 6. Relative changes (RCs) (in %) in the MPL estimates of β j under the TPS-GLM. The last two
columns indicate the AIC and R2 (Adj) of the model with dropped observations.

Dropped Obs. β0 β1 RCβ0 RCβ1 AIC R2 (Adj)

167 2.513 −0.002 0.231 0.378 1777.07 0.737
175 2.506 −0.002 0.051 1.673 1784.58 0.724
219 2.540 −0.002 1.334 7.105 1784.44 0.725
220 2.507 −0.002 0.012 0.263 1777.87 0.728
167-175 2.511 −0.002 0.169 1.052 1771.76 0.734
167-219 2.538 −0.002 1.242 6.853 1771.58 0.735
167-220 2.511 −0.002 0.179 0.884 1765.08 0.738
175-219 2.538 −0.002 1.248 7.368 1779.10 0.722
175-220 2.504 −0.002 0.104 0.684 1772.55 0.725
219-220 2.507 −0.002 0.007 0.289 1772.807 0.725
167-175-219 2.536 −0.002 1.155 7.136 1766.269 0.732
167-175-220 2.512 −0.002 0.215 3.415 1759.79 0.735
175-219-220 2.504 −0.002 0.098 0.678 1767.484 0.721
167-175-219-220 2.534 −0.002 1.072 7.800 1754.761 0.731

6. Concluding Remarks and Future Research

In this work, we study some aspects of the Thin-Plate Spline Generalized Linear Mod-
els. Specifically, we derive an iterative process to estimate the parameters and the Fisher
information matrix to approximate, through its inverse, the variance–covariance matrix
of the estimators. In addition, we extended the local influence method, obtaining closed
expressions for the Hessian and perturbation matrices under cases-weight perturbation
and additive perturbation of the response variable. We performed a statistical data analysis
with two real data sets of the agronomic and environmental area. The study showed the
advantage of incorporating a smooth surface to model the joint effect of a pair of explana-
tory variables or the spatial effect determined by the coordinates. In both applications,
Axioms 2024, 13, 346 19 of 20

it was observed that the adjusted values of the response variable were consistent. In ad-
dition, it was observed that our model presented a better fit to model the soybean yield
and ozone concentration data, compared to some classic parametric and semiparametric
models, respectively. In our analysis, it was found that those observations detected as
potentially influential generated important changes in the estimates, but not significant
inferential changes. In addition, our study confirms the need to develop the local influence
method to evaluate the sensitivity of maximum penalized likelihood estimators and thus
determine those observations that can exert an excessive influence on both the parametric
and non-parametric components, or on both.
As future work, we propose to incorporate a correlation component in the model
and extend the local influence technique to other perturbation schemes, mainly on the
non-parametric component of the model.

Author Contributions: Conceptualization, G.I.-P., P.P., M.A.U.-O. and O.N.; methodology, G.I.-P. and
O.N.; software, G.I.-P. and P.P.; validation, G.I.-P. and P.P.; formal analysis, P.P.; investigation, G.I.-P.,
P.P., O.N. and M.A.U.-O.; data curation, P.P. and M.A.U.-O.; writing—original draft preparation,
G.I.-P. and P.P.; writing—review and editing, O.N. and M.A.U.-O.; supervision, G.I.-P. and O.N.;
project administration, O.N.; funding acquisition, O.N. All authors have read and agreed to the
published version of the manuscript.
Funding: This research was funded by ANID-Fondecyt grant number 1201478.
Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.
Data Availability Statement: Data cannot be shared openly but are available on request from authors.
Acknowledgments: The third author acknowledges the support of ANID-Fondecyt 1201478.
Conflicts of Interest: The authors declare no conflicts of interest.

References
1. McCullagh, P.; Nelder, J.A. Generalized Linear Models, 2nd ed.; Chapman and Hall: London, UK, 1989.
2. Duchon, J. Interpolation des fonctions de deux variables suivant le principe de la flexion des plaques minces. RAIRO Anal. Numér.
1976, 10, 5–12. [CrossRef]
3. Duchon, J. Splines minimizing rotation-invariant semi-norms in Sobolev spaces. Lect. Notes Math. 1977, 57, 85–100.
4. Bookstein, F.L. Principal warps: Thin-plate splines and decomposition of deformations. IEEE Trans. Pattern Anal. Mach. Intell.
1989, 11, 567–585. [CrossRef]
5. Chen, C.; Li, Y.; Yan, C.; Dai, H.; Liu, G. A Thin Plate Spline-Based Feature-Preserving Method for Reducing Elevation Points
Derived from LiDAR. Remote Sens. 2015, 7, 11344–11371. [CrossRef]
6. Wahba, G. Spline Models for Observational Data; SIAM: Philadelphia, PA, USA, 1990.
7. Green, P.J.; Silverman, B.W. Nonparametric Regression and Generalized Linear Models; Chapman and Hall: Boca Raton, FL, USA, 1994.
8. Wood, S.N. Thin plate regression splines. J. R. Stat. Soc. Ser. B (Methodol.) 2003, 65, 95–114. [CrossRef]
9. Moraga, M.S.; Ibacache-Pulgar, G.; Nicolis, O. On an elliptical thin-plate spline partially varying-coefficient model. Chil. J. Stat.
2021, 12, 205–228.
10. Cook, R.D. Assessment of Local Influence. J. R. Stat. Soc. Ser. B (Methodol.) 1986, 48, 133–169. [CrossRef]
11. Thomas, W.; Cook, R.D. Assessing influence on regression coefficients in generalized linear models. Biometrika 1989, 76, 741–749.
[CrossRef]
12. Ouwens, M.N.M.; Tan, F.E.S.; Berger, M.P.F. Local influence to detect influential data structures for generalized linear mixed
models. Biometrics 2001, 57, 1166–1172. [CrossRef]
13. Zhu, H.; Lee, S. Local influence for incomplete-data models. J. R. Stat. Soc. Ser. B 2001, 63, 111–126. [CrossRef]
14. Zhu, H.; Lee, S. Local influence for generalized linear mixed models. Can. J. Stat. 2003, 31, 293–309. [CrossRef]
15. Espinheira, P.L.; Ferrari, P.L.; Cribari-Neto, F. Influence diagnostics in beta regression. Comput. Stat. Data Anal. 2008, 52, 4417–4431.
[CrossRef]
16. Rocha, A.; Simas, A. Influence diagnostics in a general class of beta regression models. TEST 2001, 20, 95–119. [CrossRef]
17. Ferrari, S.; Spinheira, P.; Cribari-Neto, F. Diagnostic tools in beta regression with varying dispersion. Stat. Neerl. 2011, 65, 337–351.
[CrossRef]
18. Ferreira, C.S.; Paula, G.A. Estimation and diagnostic for skew-normal partially linear models. J. Appl. Stat. 2017, 44, 3033–3053.
[CrossRef]
19. Emami, H. Local influence for Liu estimators in semiparametric linear models. Stat. Pap. 2018, 59, 529–544. [CrossRef]
Axioms 2024, 13, 346 20 of 20

20. Liu, Y.; Mao, G.; Leiva, V.; Liu, S.; Tapia, A. Diagnostic Analytics for an Autoregressive Model under the Skew-Normal Distribution.
Mathematics 2020, 8, 693. [CrossRef]
21. Thomas, W. Influence diagnostics for the cross-validated smoothing parameter in spline smoothing. J. Am. Stat. Assoc. 1991,
9, 693–698. [CrossRef]
22. Ibacache, G.; Paula, G.A. Local Influence for student-t partially linear models. Comput. Stat. Data Anal. 2011, 55, 1462–1478.
[CrossRef]
23. Ibacache-Pulgar, G.; Paula, G.A.; Galea, M. Influence diagnostics for elliptical semiparametric mixed models. Stat. Model. 2012,
12, 165–193. [CrossRef]
24. Ibacache, G.; Paula, G.A.; Cysneiros, F. Semiparametric additive models under symmetric distributions. Test 2013, 22, 103–121.
[CrossRef]
25. Zhang, J.; Zhang, X.; Ma, H.; Zhiya, C. Local influence analysis of varying coefficient linear model. J. Interdiscip. Math. 2015,
3, 293–306. [CrossRef]
26. Ibacache-Pulgar, G.; Reyes, S. Local influence for elliptical partially varying coefficient model. Stat. Model. 2018, 18, 149–174.
[CrossRef]
27. Ibacache-Pulgar, G.; Figueroa-Zuñiga, J.; Marchant, C. Semiparametric additive beta regression models: Inference and local
influence diagnostics. REVSTAT-Stat. J. 2019, 19, 255–274.
28. Cavieres, J.; Ibacache-Pulgar, G.; Contreras-Reyes, J. Thin plate spline model under skew-normal random errors: Estimation and
diagnostic analysis for spatial data. J. Stat. Comput. Simul. 2023, 93, 25–45. [CrossRef]
29. Jeldes, N.; Ibacache-Pulgar, G.; Marchant, C.; López-Gonzales, J.L. Modeling Air Pollution Using Partially Varying Coefficient
Models with Heavy Tails. Mathematics 2022, 10, 3677. [CrossRef]
30. Saavedra-Nievas, J.C.; Nicolis, O.; Galea, M.; Ibacache-Pulgar, G. Influence diagnostics in Gaussian spatial—Temporal linear
models with separable covariance. Environ. Ecol. Stat. 2023, 30, 131–155. [CrossRef]
31. Sánchez, L.; Ibacache-Pulgar, G.; Marchant, C.; Riquelme, M. Modeling Environmental Pollution Using Varying-Coefficients
Quantile Regression Models under Log-Symmetric Distributions. Axioms 2023, 12, 976. [CrossRef]
32. Green, P.J. Penalized Likelihood for General Semi-Parametric Regression Models. Int. Stat. Rev. 1987, 55, 245–259. [CrossRef]
33. Nelder, J.A.; Wedderburn, R.W.M. Generalized Linear Models. J. R. Stat. Soc. Ser. A (Gen.) 1972, 135, 370–384. [CrossRef]
34. Wood, S.N. Generalized Additive Models: An Introduction with R, 2nd ed.; Chapman and Hall/CRC: Boca Raton, FL, USA, 2017.
35. Akaike, H. Information theory as an extension of the maximum likelihood principle. In Proceedings of the Second International
Symposium on Information Theory; Petrov, B.N., Csaki, F., Eds.; Academiai Kiado: Budapest, Hungary, 1973.
36. Escobar, L.A.; Meeker, W.Q. Assessing Influence in Regression Analysis with Censored Data. Biometrics 1992, 48, 507–528.
[CrossRef] [PubMed]
37. Billor, N.; Loynes, R.M. Local influence: A new approach. Comm. Statist. Theory Meth. 1993, 22, 1595–1611. [CrossRef]
38. MathWorks Inc. MATLAB Version: 9.13.0 (R2022b); The MathWorks Inc.: Natick, MA, USA, 2022. Available online: https:
//www.mathworks.com (accessed on 10 October 2022).
39. Uribe-Opazo, M.A.; Borssoi, J.A.; Galea, M. Influence diagnostics in Gaussian spatial linear models. J. Appl. Stat. 2012, 3, 615–630.
[CrossRef]
40. Breiman, L.; Friedman, J.H. Estimating optimal transformations for multiple regression and correlation. J. Am. Stat. Assoc. 1985,
80, 580–598. [CrossRef]
41. Ibacache-Pulgar, G.; Lira, V.; Villegas, C. Assessing Influence on Partially Varying-coefficient Generalized Linear Model. REVSTAT-
Stat. J. 2022. Available online: https://revstat.ine.pt/index.php/REVSTAT/article/view/507 (accessed on 10 October 2022).

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual
author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to
people or property resulting from any ideas, methods, instructions or products referred to in the content.