Generalized functional linear model

The generalized functional linear model (GFLM) is an extension of the generalized linear model (GLM) that allows one to regress univariate responses of various types (continuous or discrete) on functional predictors, which are mostly random trajectories generated by a square-integrable stochastic processes. Similarly to GLM, a link function relates the expected value of the response variable to a linear predictor, which in case of GFLM is obtained by forming the scalar product of the random predictor function X {\displaystyle X} with a smooth parameter function β {\displaystyle \beta } . Functional Linear Regression, Functional Poisson Regression and Functional Binomial Regression, with the important Functional Logistic Regression included, are special cases of GFLM. Applications of GFLM include classification and discrimination of stochastic processes and functional data.^[1]

A key aspect of GFLM is estimation and inference for the smooth parameter function β {\displaystyle \beta } which is usually obtained by dimension reduction of the infinite dimensional functional predictor. A common method is to expand the predictor function X {\displaystyle X} in an orthonormal basis of L² space, the Hilbert space of square integrable functions with the simultaneous expansion of the parameter function in the same basis. This representation is then combined with a truncation step to reduce the contribution of the parameter function β {\displaystyle \beta } in the linear predictor to a finite number of regression coefficients. Functional principal component analysis (FPCA) that employs the Karhunen–Loève expansion is a common and parsimonious approach to accomplish this. Other orthogonal expansions, like Fourier expansions and B-spline expansions may also be employed for the dimension reduction step. The Akaike information criterion (AIC) can be used for selecting the number of included components. Minimization of cross-validation prediction errors is another criterion often used in classification applications. Once the dimension of the predictor process has been reduced, the simplified linear predictor allows to use GLM and quasi-likelihood estimation techniques to obtain estimates of the finite dimensional regression coefficients which in turn provide an estimate of the parameter function β {\displaystyle \beta } in the GFLM.

The predictor functions X ( t ) , t ∈ T {\displaystyle \textstyle X(t),t\in T} , typically are square integrable stochastic processes on a real interval T {\displaystyle T} and the unknown smooth parameter function β ( t ) , t ∈ T {\displaystyle \beta (t),t\in T} , is assumed to be square integrable on T {\displaystyle T} . Given a real measure d w {\displaystyle dw} on T {\displaystyle T} , the linear predictor is given by η = α + ∫ X c ( t ) β ( t ) d w ( t ) {\displaystyle \eta =\alpha +\int X^{c}(t)\beta (t)\,dw(t)} where X c ( t ) = X ( t ) − E ( X ( t ) ) {\displaystyle X^{c}(t)=X(t)-{\text{E}}(X(t))} is the centered predictor process and α {\displaystyle \alpha } is a scalar that serves as an intercept.

The outcome Y {\displaystyle Y} is typically a real valued random variable which may be either continuous or discrete. Often the conditional distribution of Y {\displaystyle Y} given the predictor process is specified within the exponential family. However it is also sufficient to consider the functional quasi-likelihood set up, where instead of the distribution of the response one specifies the conditional variance function, V a r ( Y ∣ X ) = σ 2 ( μ ) {\displaystyle {\rm {{Var}(Y\mid X)=\sigma ^{2}(\mu )}}} , as a function of the conditional mean, E ( Y ∣ X ) = μ {\displaystyle {\rm {{E}(Y\mid X)=\mu }}} .

The link function g {\displaystyle g} is a smooth invertible function, that relates the conditional mean of the response E ( Y ∣ X ) = μ {\displaystyle {\rm {{E}(Y\mid X)=\mu }}} with the linear predictor η = α + ∫ X c ( t ) β ( t ) d w ( t ) {\displaystyle \eta =\alpha +\int X^{c}(t)\beta (t)\,dw(t)} . The relationship is given by μ = g ( η ) {\displaystyle \mu =g(\eta )} .

In order to implement the necessary dimension reduction, the centered predictor process X c ( t ) {\displaystyle X^{c}(t)} and the parameter function β ( t ) {\displaystyle \beta (t)} are expanded as,

X c ( t ) = ∑ j = 1 ∞ ξ j ρ j ( t )  and  β ( t ) = ∑ j = 1 ∞ β j ρ j ( t ) , {\displaystyle X^{c}(t)=\sum _{j=1}^{\infty }\xi _{j}\rho _{j}(t){\text{ and }}\beta (t)=\sum _{j=1}^{\infty }\beta _{j}\rho _{j}(t),}

where ρ j , j = 1 , 2 , … {\displaystyle \rho _{j},j=1,2,\ldots } is an orthonormal basis of the function space L 2 ( d w ) , {\displaystyle L^{2}(dw),} such that ∫ T ρ j ( t ) ρ k ( t ) d w ( t ) = δ j k {\displaystyle \int _{T}\rho _{j}(t)\rho _{k}(t)\,dw(t)=\delta _{jk}} where δ j k = 1 {\displaystyle \delta _{jk}=1} if j = k {\displaystyle j=k} and 0 {\displaystyle 0} otherwise.

The random variables ξ j {\displaystyle \xi _{j}} are given by ξ j = ∫ X c ( t ) ρ j ( t ) d w ( t ) {\displaystyle \xi _{j}=\int X^{c}(t)\rho _{j}(t)\,dw(t)} and the coefficients β j {\displaystyle \beta _{j}} as β j = ∫ β ( t ) ρ j ( t ) d w ( t ) {\displaystyle \beta _{j}=\int \beta (t)\rho _{j}(t)\,dw(t)} for j = 1 , 2 , … {\displaystyle j=1,2,\ldots } .

E ( ξ j ) = 0 {\displaystyle {\text{E}}(\xi _{j})=0} and ∑ j = 1 ∞ β j 2 < ∞ {\displaystyle \sum _{j=1}^{\infty }\beta _{j}^{2}<\infty } and denoting σ j 2 = Var ( ξ j ) = E ( ξ j 2 ) {\displaystyle \sigma _{j}^{2}={\text{Var}}(\xi _{j})={\text{E}}(\xi _{j}^{2})} , so ∑ j = 1 ∞ σ j 2 = ∫ E ( X c ( t ) ) 2 d w ( t ) < ∞ {\displaystyle \sum _{j=1}^{\infty }\sigma _{j}^{2}=\int {\text{E}}(X^{c}(t))^{2}\,dw(t)<\infty } .

From the orthonormality of the basis functions ρ j {\displaystyle \rho _{j}} , it follows immediately that ∫ X c ( t ) β ( t ) d w ( t ) = ∑ j = 1 ∞ β j ξ j {\displaystyle \int X^{c}(t)\beta (t)\,dw(t)=\sum _{j=1}^{\infty }\beta _{j}\xi _{j}} .

The key step is then approximating η = α + ∫ X c ( t ) β ( t ) d w ( t ) = α + ∑ j = 1 ∞ β j ξ j {\displaystyle \eta =\alpha +\int X^{c}(t)\beta (t)\,dw(t)=\alpha +\sum _{j=1}^{\infty }\beta _{j}\xi _{j}} by η ≈ α + ∑ j = 1 p β j ξ j {\displaystyle \eta \approx \alpha +\sum _{j=1}^{p}\beta _{j}\xi _{j}} for a suitably chosen truncation point p {\displaystyle p} .

FPCA gives the most parsimonious approximation of the linear predictor for a given number of basis functions as the eigenfunction basis explains more of the variation than any other set of basis functions.

For a differentiable link function with bounded first derivative, the approximation error of the p {\displaystyle p} -truncated model i.e. the linear predictor truncated to the summation of the first p {\displaystyle p} components, is a constant multiple of Var ( ∑ j = p + 1 ∞ β j ξ j ) = E ( ( ∑ j = p + 1 ∞ β j ξ j ) 2 ) = ∑ j = p + 1 ∞ β j σ j 2 {\displaystyle {\text{Var}}(\sum _{j=p+1}^{\infty }\beta _{j}\xi _{j})={\text{E}}\left(\left(\sum _{j=p+1}^{\infty }\beta _{j}\xi _{j}\right)^{2}\right)=\sum _{j=p+1}^{\infty }\beta _{j}\sigma _{j}^{2}} .

A heuristic motivation for the truncation strategy derives from the fact that E ( ( ∑ j = p + 1 ∞ β j ξ j ) 2 ) = ∑ j = p + 1 ∞ β j σ j 2 ≤ ∑ j = p + 1 ∞ β j 2   ∑ j = p + 1 ∞ σ j 2 {\displaystyle {\text{E}}\left(\left(\sum _{j=p+1}^{\infty }\beta _{j}\xi _{j}\right)^{2}\right)=\sum _{j=p+1}^{\infty }\beta _{j}\sigma _{j}^{2}\leq \sum _{j=p+1}^{\infty }\beta _{j}^{2}\ \sum _{j=p+1}^{\infty }\sigma _{j}^{2}} which is a consequence of the Cauchy–Schwarz inequality and by noting that the right hand side of the last inequality converges to 0 as p → ∞ {\displaystyle p\rightarrow \infty } since both ∑ j = 1 ∞ β j 2 {\displaystyle \sum _{j=1}^{\infty }\beta _{j}^{2}} and ∑ j = 1 ∞ σ j 2 {\displaystyle \sum _{j=1}^{\infty }\sigma _{j}^{2}} are finite.

For the special case of the eigenfunction basis, the sequence σ j 2 , j = 1 , 2 , … {\displaystyle \sigma _{j}^{2},j=1,2,\ldots } corresponds to the sequence of the eigenvalues of the covariance kernel G ( s , t ) = Cov ( X ( s ) , X ( t ) ) ,   s , t ∈ T {\displaystyle G(s,t)={\text{Cov}}(X(s),X(t)),\ s,t\in T} .

For data with n {\displaystyle n} i.i.d observations, setting ξ j 0 = 1 {\displaystyle \xi _{j}^{0}=1} , β 0 = α {\displaystyle \beta _{0}=\alpha } and ξ j i = ∫ X i ( t ) ρ j ( t ) d w ( t ) {\displaystyle \xi _{j}^{i}=\int X_{i}(t)\rho _{j}(t)\,dw(t)} , the approximated linear predictors can be represented as η i = ∑ j = 0 p β j ξ j i , i = 1 , 2 , … , n {\displaystyle \eta _{i}=\sum _{j=0}^{p}\beta _{j}\xi _{j}^{i},i=1,2,\ldots ,n} which are related to the means through μ i = g ( η i ) {\displaystyle \mu _{i}=g(\eta _{i})} .

The main aim is to estimate the parameter function β {\displaystyle \beta } .

Once p {\displaystyle p} has been fixed, standard GLM and quasi-likelihood methods can be used for the p {\displaystyle p} -truncated model to estimate β T = ( β 0 , β 1 , … , β p ) {\displaystyle {\boldsymbol {\beta }}^{T}=(\beta _{0},\beta _{1},\ldots ,\beta _{p})} by solving the estimating equation or the score equation U ( β ) = 0. {\displaystyle U(\beta )=0.}

The vector valued score function turns out to be U ( β ) = ∑ i = 1 n ( Y i − μ i ) g ′ ( η i ) ξ i / σ 2 ( μ i ) {\displaystyle U(\beta )=\sum _{i=1}^{n}(Y_{i}-\mu _{i})g'(\eta _{i})\xi _{i}/\sigma ^{2}(\mu _{i})} which depends on β {\displaystyle {\boldsymbol {\beta }}} through μ {\displaystyle \mu } and η {\displaystyle \eta } .

Just as in GLM, the equation U ( β ) = 0 {\displaystyle U(\beta )=0} is solved using iterative methods like Newton–Raphson (NR) or Fisher scoring (FS) or iteratively reweighted least squares (IWLS) to get the estimate of the regression coefficients β ^ {\displaystyle {\boldsymbol {\hat {\beta }}}} , leading to the estimate of the parameter function β ^ ( t ) = β ^ o + ∑ j = 1 p β ^ j ρ j ( t ) {\displaystyle {\hat {\beta }}(t)={\hat {\beta }}_{o}+\sum _{j=1}^{p}{\hat {\beta }}_{j}\rho _{j}(t)} . When using the canonical link function, these methods are equivalent.

Results are available in the literature of p {\displaystyle p} -truncated models as p → ∞ {\displaystyle p\rightarrow \infty } which provide asymptotic inference for the deviation of the estimated parametric function from the true parametric function and also asymptotic tests for regression effects and asymptotic confidence regions.

If the response variable Y i {\displaystyle Y_{i}} , given X i ∈ L 2 ( T ) {\displaystyle X_{i}\in L^{2}(T)} follows the one parameter exponential family, then its probability density function or probability mass function (as the case may be) is

f ( y i ∣ X i ) = exp ⁡ ( y i θ i − b ( θ i ) ϕ + c ( y i , ϕ ) ) {\displaystyle f(y_{i}\mid X_{i})=\exp \left({\frac {y_{i}\theta _{i}-b(\theta _{i})}{\phi }}+c(y_{i},\phi )\right)}

for some functions b {\displaystyle b} and c {\displaystyle c} , where θ i {\displaystyle \theta _{i}} is the canonical parameter, and ϕ {\displaystyle \phi } is a dispersion parameter which is typically assumed to be positive.

In the canonical set up, η i = α + ∫ X i c ( t ) β ( t ) d w ( t ) = θ i {\displaystyle \eta _{i}=\alpha +\int X_{i}^{c}(t)\beta (t)\,dw(t)=\theta _{i}} and from the properties of exponential family,

μ i = b ′ ( θ i ) ,  and so  μ i = b ′ ( η i ) . {\displaystyle \mu _{i}=b'(\theta _{i}),{\text{ and so }}\mu _{i}=b'(\eta _{i}).}

Hence b ′ {\displaystyle b'} serves as a link function and is called the canonical link function.

Var ( y i ) = ϕ b ″ ( θ i ) = ϕ b ″ ( η i ) = ϕ g ′ ( η i ) = ϕ g ′ ( g − 1 ( μ i ) ) ) {\displaystyle {\text{Var}}(y_{i})=\phi b''(\theta _{i})=\phi b''(\eta _{i})=\phi g'(\eta _{i})=\phi g'(g^{-1}(\mu _{i})))} is the corresponding variance function and ϕ {\displaystyle \phi } the dispersion parameter.

Functional linear regression, one of the most useful tools of functional data analysis, is an example of GFLM where the response variable is continuous and is often assumed to have a Normal distribution. The variance function is a constant function and the link function is identity. Under these assumptions the GFLM reduces to the FLR,

μ = E ⁡ ( Y ∣ X ) = η = α + ∫ X c ( t ) β ( t ) d w ( t ) {\displaystyle \mu =\operatorname {E} (Y\mid X)=\eta =\alpha +\int X^{c}(t)\beta (t)\,dw(t)}

Without the normality assumption, the constant variance function motivates the use of quasi-normal techniques.

When the response variable has binary outcomes, i.e., 0 or 1, the distribution is usually chosen as Bernoulli, and then μ i = P ( Y i = 1 ∣ X i ) {\displaystyle \mu _{i}=P(Y_{i}=1\mid X_{i})} . Popular link functions are the expit function, which is the inverse of the logit function (functional logistic regression) and the probit function (functional probit regression). Any cumulative distribution function F has range [0,1] which is the range of binomial mean and so can be chosen as a link function. Another link function in this context is the complementary log–log function, which is an asymmetric link. The variance function for binary data is given by Var ⁡ ( Y i ) = ϕ μ i ( 1 − μ i ) {\displaystyle \operatorname {Var} (Y_{i})=\phi \mu _{i}(1-\mu _{i})} where the dispersion parameter ϕ {\displaystyle \phi } is taken as 1 or alternatively the quasi-likelihood approach is used.

Another special case of GFLM occurs when the outcomes are counts, so that the distribution of the responses is assumed to be Poisson. The mean μ i {\displaystyle \mu _{i}} is typically linked to the linear predictor η i {\displaystyle \eta _{i}} via a log-link, which is also the canonical link . The variance function is Var ⁡ ( Y i ) = ϕ μ i {\displaystyle \operatorname {Var} (Y_{i})=\phi \mu _{i}} , where the dispersion parameter ϕ {\displaystyle \phi } is 1, except when the data might be over-dispersed which is when the quasi-Poisson approach is used.

Extensions of GFLM have been proposed for the cases where there are multiple predictor functions.^[2] Another generalization is called the Semi Parametric Quasi-likelihood Regression (SPQR)^[1] which considers the situation where the link and the variance functions are unknown and are estimated non-parametrically from the data. This situation can also be handled by single or multiple index models, using for example Sliced Inverse Regression (SIR).

Another extension in this domain is Functional Generalized Additive Model (FGAM))^[3] which is a generalization of generalized additive model(GAM) where

g − 1 ( E ⁡ ( Y ∣ X ) ) = α + ∑ j = 1 p f j ( ξ j ) , {\displaystyle g^{-1}(\operatorname {E} (Y\mid X))=\alpha +\sum _{j=1}^{p}f_{j}(\xi _{j}),}

where ξ j {\displaystyle \xi _{j}} are the expansion coefficients of the random predictor function X {\displaystyle X} and each f j {\displaystyle f_{j}} is an unknown smooth function that has to be estimated and where E ( f j ( ξ j ) ) = 0. {\displaystyle {\text{E}}(f_{j}(\xi _{j}))=0.} .

In general, estimation in FGAM requires combining IWLS with backfitting. However, if the expansion coefficients are obtained as functional principal components, then in some cases (e.g. Gaussian predictor function X {\displaystyle X} ), they will be independent in which case backfitting is not needed, and one can use popular smoothing methods for estimating the unknown parameter functions f j {\displaystyle f_{j}} .

A popular data set that has been used for a number of analysis in the domain of functional data analysis consists of the number of eggs laid daily until death of 1000 Mediterranean fruit flies (or medflies for short)[1][2]. The plot here shows the egg laying trajectories in the first 25 days of life of about 600 female medflies (those that have at least 20 remaining eggs in their lifetime). The red colored curves belong to those flies that will lay less than the median number of remaining eggs, while the blue colored curves belong to the flies that will lay more than the median number of remaining eggs after age 25. A related problem of classifying medflies as long-lived or short-lived based on the initial egg laying trajectories as predictors and the subsequent longevity of the flies as response has been studied with the GFLM^[1]