Manuel Oviedo de la Fuente

Functional linear model (FLR) with basis representation

Supose that \(\mathcal{X} \in \mathcal{L}_{2}(T)\) and [[“y”]]$. Assume also that $[(t)]=0, t and [y]=0.

The FLM states that \[y= \left\langle \mathcal{X},\beta \right\rangle +\varepsilon=\int_{T}X(t)\beta(t)dt+\varepsilon\] where \(\beta \in \mathcal{L}_{2}(T)\) and \(\varepsilon\) is the errror term.

One way of estimating \(\beta\), it is representing the parametmer (and \(\mathcal{X}\)) in a \(\mathcal{L}_2\)-basis in the following way:

\[\beta(t)=\sum_k \beta_k \theta_k(t), \mathbf{X}(t)=\sum_k c_i \psi_k(t)\]

• fregre.basis() fucntion uses fixed basis: B–spline, Fourier, etc. Ramsay and Silverman (2005), Cardot, Ferraty, and Sarda (1999))

The next code illustrates how to estimate the fat contents using a sample of absorbances curves.

rangett <- X$rangeval
basis1 = create.bspline.basis(rangeval = rangett, nbasis = 17)
basis2 = create.bspline.basis(rangeval = rangett, nbasis = 7)
res.basis0 = fregre.basis(X, y, basis.x = basis1, basis.b = basis2)
res.basis1 = fregre.basis(X.d1, y, basis.x = basis1, basis.b = basis2)
res.basis2 = fregre.basis(X.d2, y, basis.x = basis1, basis.b = basis2)
res.basis0$r2;res.basis1$r2;res.basis2$r2

## [1] 0.9521075

## [1] 0.9484766

## [1] 0.9583489

summary(res.basis2)

##  *** Summary Functional Data Regression with representation in Basis *** 
## 
## Call:
## fregre.basis(fdataobj = X.d2, y = y, basis.x = basis1, basis.b = basis2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.3174  -1.6883  -0.1376   1.7609   5.7993 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.750e+01  2.416e-01  72.417  < 2e-16 ***
## X.d2.bspl4.1 -1.732e+04  3.396e+03  -5.098 1.28e-06 ***
## X.d2.bspl4.2  1.316e+04  2.658e+03   4.951 2.42e-06 ***
## X.d2.bspl4.3 -1.383e+03  1.320e+03  -1.047   0.2970    
## X.d2.bspl4.4  2.478e+02  9.591e+02   0.258   0.7966    
## X.d2.bspl4.5 -2.707e+02  1.000e+03  -0.271   0.7872    
## X.d2.bspl4.6  3.346e+03  1.696e+03   1.973   0.0508 .  
## X.d2.bspl4.7 -5.329e+03  1.298e+03  -4.105 7.38e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.744 on 121 degrees of freedom
## Multiple R-squared:  0.9583, Adjusted R-squared:  0.9559 
## F-statistic: 397.7 on 7 and 121 DF,  p-value: < 2.2e-16
## 
## -Names of possible atypical curves: 43 
## -Names of possible influence curves: 34 140 45

par(mfrow=c(1,3))
plot(res.basis0$beta.est)

## [1] "done"

plot(res.basis1$beta.est)

## [1] "done"

plot(res.basis2$beta.est)

## [1] "done"

The choice of the appropiate basis (and the number of basis elements) becomes now in a crucial step:

res.basis.cv = fregre.basis(X, y)
summary(res.basis.cv)

##  *** Summary Functional Data Regression with representation in Basis *** 
## 
## Call:
## fregre.basis(fdataobj = X, y = y)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.1454  -1.6915  -0.0624   1.8437   5.5219 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   17.4969     0.2378  73.571  < 2e-16 ***
## X.bspl4.1   -211.7506    66.2778  -3.195 0.001795 ** 
## X.bspl4.2    275.8737    78.5113   3.514 0.000627 ***
## X.bspl4.3   -175.4673    50.2726  -3.490 0.000679 ***
## X.bspl4.4     54.1873    23.6024   2.296 0.023451 *  
## X.bspl4.5     -7.9114    13.8454  -0.571 0.568809    
## X.bspl4.6     -2.0064    11.6901  -0.172 0.864017    
## X.bspl4.7     11.3758    19.6239   0.580 0.563228    
## X.bspl4.8    -33.8397    48.2034  -0.702 0.484050    
## X.bspl4.9     54.1517    81.8808   0.661 0.509679    
## X.bspl4.10   -42.9837    72.3735  -0.594 0.553706    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.701 on 118 degrees of freedom
## Multiple R-squared:  0.9606, Adjusted R-squared:  0.9573 
## F-statistic:   288 on 10 and 118 DF,  p-value: < 2.2e-16
## 
## -Names of possible atypical curves: 43 
## -Names of possible influence curves: 34 140 6 7 45 86

• Functional Principal Components (FPC).(Cardot, Ferraty, and Sarda 1999), fregre.pc()

x<-X
basis.pc0 = create.pc.basis(X,1:3)
res.pc1 = fregre.pc(X, y, basis.x = basis.pc)
summary(res.pc1)

##  *** Summary Functional Data Regression with Principal Components ***
## 
## Call:
## fregre.pc(fdataobj = X, y = y, basis.x = basis.pc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -24.6464  -4.4747  -0.1008   6.1505  15.6592 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  17.4969     0.7620  22.962  < 2e-16 ***
## PC1          -0.9465     0.1074  -8.814 8.65e-15 ***
## PC2           2.5302     1.0578   2.392   0.0183 *  
## PC3          18.7628     1.8082  10.377  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.654 on 125 degrees of freedom
## Multiple R-squared:  0.572,  Adjusted R-squared:  0.5618 
## F-statistic:  55.7 on 3 and 125 DF,  p-value: < 2.2e-16
## 
## 
## -With 3 Principal Components is  explained  99.45 %
##  of the variability of explicative variables. 
## 
## -Variability for each  principal components -PC- (%):
##   PC1   PC2   PC3 
## 98.12  0.99  0.35 
## -Names of possible atypical curves: No atypical curves 
## -Names of possible influence curves: 34 140 99 139

res.pc2 = fregre.pc.cv(X, y)
summary(res.pc2)

##            Length Class     Mode   
## fregre.pc  19     fregre.fd list   
## pc.opt      6     -none-    numeric
## lambda.opt  1     -none-    numeric
## PC.order    8     -none-    numeric
## MSC.order   8     -none-    numeric

par(mfrow=c(1,2))
plot(res.pc1$beta.est)
plot(res.pc2[[1]]$beta.est)

FLM with functional and non functional covariates

\[E(y)=\alpha+\mathbf{Z}\beta+\sum_{q=1}^Q \left\langle \mathcal{X}^{q}(t),\beta_{q}(t)\right\rangle \]

where \(\left\{\mathcal{X}_q(t)\right\}_{q=1}^Q\) are function covariates and \(\mathbf{Z}=\left\{{Z_j}\right\}_{j=1}^J\) the non–functional covariates.

dataf = as.data.frame(tecator[["y"]][ind,]) # Fat, Protein, Water
basis.pc2 = create.pc.basis(X.d2,1:4)
basis.x = list(X = basis.pc0, X.d2 =basis.pc2)
f = Fat ~ X+X.d2
ldata = list(df = dataf, X=X,X.d2=X.d2)
res.lm1 = fregre.lm(f, ldata, basis.x = basis.x)
f = Fat ~ Water+X.d2
res.lm2 = fregre.lm(f, ldata, basis.x = basis.x)

## 
## Call:
## lm(formula = pf, data = XX, x = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.7919  -1.6743  -0.0009   1.9080   7.0892 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.750e+01  2.519e-01  69.468  < 2e-16 ***
## X.PC1       -2.997e-02  7.146e-02  -0.419  0.67570    
## X.PC2       -1.246e+00  2.356e+00  -0.529  0.59791    
## X.PC3        9.634e+00  5.301e+00   1.817  0.07164 .  
## X.d2.PC1     3.962e+03  2.788e+02  14.213  < 2e-16 ***
## X.d2.PC2    -7.007e+02  2.128e+03  -0.329  0.74255    
## X.d2.PC3     9.204e+02  5.999e+02   1.534  0.12755    
## X.d2.PC4    -7.759e+03  2.333e+03  -3.326  0.00117 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.861 on 121 degrees of freedom
## Multiple R-squared:  0.9547, Adjusted R-squared:  0.9521 
## F-statistic: 364.6 on 7 and 121 DF,  p-value: < 2.2e-16

## [1] "No newx entered"

Other procedures

4. Other procedures

• Partial Least Squares (FPLS). fregre.pls(), (Krämer and Sugiyama 2011, @pedra2005)
• Penalized versions and parameter selection: fregre.pc.cv, fregre.basis.cv, fregre.np.cv (Febrero-Bande and Oviedo de la Fuente 2012)
• F-test for the FLM with scalar response: flm.Ftest, F-test (Garcı́a-Portugués, González-Manteiga, and Febrero-Bande 2014)
• Goodness-of-fit test for the FLM with scalar response: flm.test (Garcı́a-Portugués, González-Manteiga, and Febrero-Bande 2014)
• Measures of influence in FLM with scalar response: influence.fdata,(Febrero-Bande, Galeano, and González-Manteiga 2010)
• Beta parameter estimation by wild or smoothed bootstrap procedure: fregre.bootstrap
• FLM with a functional response: fregre.basis.fr (Chiou et al. 2004)

Back to Top

Non Linear Model (Ferraty and Vieu 2006)

Supose \((\mathcal{X},Y)\) are a pair of r.v. with \(y\in \mathbb{R}\) where \(\mathbb{E}\) is a semi-metric space. To predict the resonse \(Y\) with \(\mathcal{X}\), the estimation is:

\[m(\mathcal{X})=\mathbb{E}(Y|X=\mathcal{X})\], where the NW estimator is given by:

\[\hat{m}(\mathcal{X})=\frac{\sum_{i=1}^n Y_i{ K(d(\mathcal{X},X_i)/h)}}{\sum_{i=1}^n {K(d(\mathcal{X},X_i)/h)}}\]

where K is an asymmetric kernel function and h is the bandwidth parameter.

Semi Linear Model (Aneiros-Pérez and Vieu 2006)

Let \((\mathcal{X},\mathbf{Z},y)\) with \(y\in \mathbb{R}\) (response), \(\mathcal{X}\in \mathbb{E}\) (functional) and \(\mathbf{Z} \in \mathbb{R}^p\) (MV covariates).

\[y = Z + m(X) + \varepsilon\]

Arguments for fregre.np() and fregre.plm() function

• –Ker–: type of asymmetric kernel function, by default asymmetric normal kernel (cosine, epanechnicov, quadratic,….).
• –metric–: type of metric or semimetric. –type.S–: type of smoothing matrix \(\mathbf{S}\): S.NW, S.LLR, S.KNN.

tecator<-list("df"=tecator$y,"absorp.fdata"=tecator$absorp.fdata)
X=tecator$absorp.fdata
y<-tecator$df$Fat

np<-fregre.np(X, y, metric = semimetric.deriv, nderiv = 1,type.S = S.KNN)

Again, it has also implemented the function fregre.np.cv to estimate the smoothing parameter \(h\) by the validation criteria.

np<-fregre.np(X, y, metric = semimetric.deriv, nderiv = 1,type.S = S.KNN)
np.cv<-fregre.np.cv(X, y, metric = semimetric.deriv, nderiv = 1,type.S = S.KNN,h=c(3:9))
c(np$h.opt,np.cv$h.opt)

## [1] 6 3

c(np$r2,np.cv$r2)

## [1] 0.8887073 0.9515189

Back to Top

Generalized Linear Models (Müller and Stadtmüller 2005)

One natural extension of LM model is the generalized functional linear regression model (GFLM) which allows various types of the response. In the GLM framework it is generally assumed that \(y_i|X_i\) can be chosen within the set of distributions belonging to the exponential family.

In Generalized Functional Linear Model (FGLM), The scalar response \(y\)(belonging to a Exponential Family PDF) is estimated by functional \(\left\{\mathcal{X}_q(t)\right\}_{q=1}^Q\) and also non–functional \(\mathbf{Z}=\left\{{Z_j}\right\}_{j=1}^J\) covariates by:

\[E(y)=g^{-1}\left(\alpha+\mathbf{Z}\beta+\sum_{q=1}^Q \left\langle \mathcal{X}^{q}(t),\beta_{q}(t)\right\rangle\right) \] where \(g()\) is the inverse link function.

Example of logistic regression

In logistic regression, the probability, \(\pi_i\) , of the occurrence of an event, \(Y_i = 1\), rather than the event \(Y_i = 0\), conditional on a vector of covariates \(\mathcal{X}_i(t)\) is expressed as:

\[ p_i = \mathbb{P}[Y = 1|{X_i(t): t \in T }]=\frac{+exp\left\{\alpha+\int_{T}X_{i}(t)\beta(t)dt \right\}}{1+exp\left\{\alpha+\int_{T}X_{i}(t)\beta(t)dt \right\}}\ , i= 1,\ldots,n\]

with \(\epsilon\) are the independent errors with zero mean.

data(tecator)
names(tecator)[2]<-"df"
tecator$df$fat15<-ifelse(tecator$df$Fat<15,0,1)
tecator$absorp.d2=fdata.deriv(tecator$absorp.fdata,nderiv=2)
res.glm<-fregre.glm(fat15 ~ absorp.d2,data=tecator,family=binomial())
#summary(a)
yfit<-ifelse(res.glm$fitted.values<.5,0,1)
table(tecator$df$fat15,yfit)

##    yfit
##       0   1
##   0 109   3
##   1   2 101

Back to Top

Generalized Functional Additive Model

1. Generalized Functional Spectral Additive Linear Model (FGSAM), (Müller and Yao 2012)

\[E(y)=g^{-1}\left(\alpha+\sum_{j=1}^J f_{j}\left(\mathbf{Z}^{j}\right)+\sum_{q=1}^Q s_q\left(\mathcal{X}_{i}^{q}(t)\right)\right)\]

where \({f}(\cdot),{s}(\cdot)\) are the smoothed functions.

res.gsam<-fregre.gsam(fat15~ s(absorp.d2),data=tecator,family=binomial())
yfit<-ifelse(res.gsam$fitted<.5,0,1)
table(tecator$df$fat15,yfit)

##    yfit
##       0   1
##   0 112   0
##   1   0 103

2. Generalized Functional Kernel Additive Linear Model (FGKAM), (Febrero-Bande and González-Manteiga 2013)

\[E(y)=g^{-1}\left(\alpha+\sum_{q=1}^Q\mathcal{K}\left(\mathcal{X}^{q}_i(t)\right)\right)\] where \(\mathcal{K}(\cdot)\) is the kernel estimator.

tecator2<-tecator[-1]
res.gkam<-fregre.gkam(fat15 ~ absorp.d2,data=tecator2,family=binomial(),
             control = list(maxit = 10))
res.gkam

## 
## Family: binomial 
## Link function: logit 
## 
## alpha= -2.47   n=  215 
## ****    ****    ****    ****    ****    ****
##                     h cor(f(X),eta)  edf
## f(absorp.d2) 0.000469         0.957 64.7
## ****    ****    ****    ****    ****    ****
## edf: Equivalent degrees of freedom
## AIC= 148.354   Deviance explained = 94.3 % 
## R-sq.= 0.963   R-sq.(adj)= 0.948

yfit<-ifelse(res.gkam$fitted.values<.5,0,1)
table(tecator$df$fat15,yfit)

##    yfit
##       0   1
##   0 112   0
##   1   2 101

Back to Top

Functional GLS model

See Oviedo de la Fuente et al. (2018) for more details about the below algorithm:

A. Jointly estimation (nlme package): Minimize for \((\beta,\theta)\) the GLS criteria, i.e,

\[\Psi(\beta,\theta)=\left(y-\left\langle X,\beta \right\rangle\right)\Sigma(\theta)^{-1}\left(y-\left\langle X,\beta \right\rangle\right)\]

B. Iterative Estimation: In multivariate case, Zivot and Wang (2007) show that estimation of \(\beta\) by \(\hat{\beta}_{ML}\) is equivalent to the iterative estimation of \(\hat{\beta}\) recomputed at each iteration by the update estimator of \(\Sigma\).

1. Begin with a preliminary estimation of \(\hat{\theta}=\theta_0\). Compute \(\hat{W}=\Sigma(\theta_0)^{-1}\).

2. Estimate \({b}_\Sigma={(Z^\prime\hat{W}Z)^{-1}Z^\prime\hat{W}}y\)

3. Based on \(\hat{e}=({y-{Z}{b}_\Sigma})\), update \(\hat{\theta}=\rho({\hat{e}})\) where \(\rho\) depends on the dependence structure chosen.

4. Repeat steps 2 and 3 until convergence.

The generalized correlated cross-validation (GCCV) criterion is an extension to GCV within the context of correlated errors, Carmack, Spence, and Schucany (2012). It is defined as follows:

\[GCCV(K_x,K_\beta,\mathbf{b},\phi)=\frac{\sum_{i=1}^n \left(y_{i}-\hat{y}_{i,\mathbf{b}}\right)^2}{ \left({1-\frac{{tr}(\mathbf{G})}{n}}\right)^2} \]

where \({G}=2{H}\Sigma(\phi)-{H}\Sigma(\phi)H^\prime\) takes into account the effect of the dependence, the trace of \({G}\) is an estimation of the degrees of freedom consumed by the model and \({H}\) is the hat matrix.

The important advantage of this criterion is that it is rather easy to compute because it avoids the need to compute the inverse of the matrix \(\Sigma\). Even so, the complexity of the GLS criterion depends on the structure of \(\Sigma\) and it could sometimes be hard either to minimize or computationally expensive.

data(tecator)
ts.plot(tecator[["y"]][,"Fat"])

cor(tecator[["y"]][,"Fat",drop=F],tecator[["y"]][,"Water",drop=F])

##          Water
## Fat -0.9881002

cor(tecator[["y"]][,"Fat",drop=F],tecator[["y"]][,"Protein",drop=F])

##        Protein
## Fat -0.8608965

dcor.xy(tecator[["y"]][,"Fat",drop=F],tecator[["y"]][,"Water",drop=F])

## 
##  dcor t-test of independence
## 
## data:  D1 and D2
## T = 571.71, df = 22789, p-value < 2.2e-16
## sample estimates:
## Bias corrected dcor 
##           0.9668619

dcor.xy(tecator[["y"]][,"Fat",drop=F],tecator[["y"]][,"Protein",drop=F])

## 
##  dcor t-test of independence
## 
## data:  D1 and D2
## T = 155.43, df = 22789, p-value < 2.2e-16
## sample estimates:
## Bias corrected dcor 
##           0.7173448

x.d2<-fdata.deriv(tecator[["absorp.fdata"]],nderiv=2)
ldata=list("df"=tecator[["y"]],"x.d2"=x.d2)

res.gls=fregre.gls(Fat~x.d2, data=ldata, correlation=corAR1())
coef(res.gls[["modelStruct"]],F)

## corStruct.Phi 
##     0.4942661

res.igls=fregre.igls(Fat~x.d2, data=ldata, correlation=list("cor.ARMA"=list()),control=list("p"=1))
coef(res.igls[["corStruct"]][[1]])

##       ar1 
## 0.4889102

res.igls

## 
## Call:
## "fregre.basis"
## 
## Coefficients:
##  (Intercept)  x.d2.bspl4.1  x.d2.bspl4.2  x.d2.bspl4.3  x.d2.bspl4.4  
##        18.12       -608.07       6202.92      -8252.57       6271.33  
## x.d2.bspl4.5  
##     -7156.75

res.igls$corStruct

## $ar
## 
## Call:
## arima(x = x, order = c(p, d, q), include.mean = FALSE, transform.pars = TRUE)
## 
## Coefficients:
##          ar1
##       0.4889
## s.e.  0.0600
## 
## sigma^2 estimated as 8.076:  log likelihood = -529.76,  aic = 1063.53

data(aemet)
log10precfdata<-aemet$logprec; tempfdata<-aemet$temp
res2<-fregre.basis.fr(tempfdata,log10precfdata)
i<-1

plot(log10precfdata[i],lty=1,main=paste0("Weather station, ",i))
lines(res2$fitted.values[i],lty=2,lwd=2,col=4)