Type:	Package
Title:	Functional Data Analysis and Utilities for Statistical Computing
Version:	2.2.0
Date:	2024-11-04
Depends:	R (≥ 3.5.0), fda, splines, MASS, mgcv, knitr
Imports:	methods, grDevices, graphics, utils, stats, nlme, doParallel, parallel, iterators, foreach, kSamples
Suggests:	rmarkdown
Description:	Routines for exploratory and descriptive analysis of functional data such as depth measurements, atypical curves detection, regression models, supervised classification, unsupervised classification and functional analysis of variance.
License:	GPL-2
URL:	https://github.com/moviedo5/fda.usc, https://moviedo5.github.io/fda.usc/, https://www.jstatsoft.org/v51/i04/
BugReports:	https://github.com/moviedo5/fda.usc/issues
LazyLoad:	yes
Author:	Manuel Febrero Bande [aut], Manuel Oviedo de la Fuente [aut, cre], Pedro Galeano [ctb], Alicia Nieto [ctb], Eduardo Garcia-Portugues [ctb]
Maintainer:	Manuel Oviedo de la Fuente <manuel.oviedo@udc.es>
NeedsCompilation:	yes
Repository:	CRAN
Encoding:	UTF-8
RoxygenNote:	7.3.2
Packaged:	2024-11-09 18:10:53 UTC; Manuel Oviedo
Date/Publication:	2024-11-09 18:40:02 UTC

Functional Data Analysis and Utilities for Statistical Computing (fda.usc)

Description

This devel version carries out exploratory and descriptive analysis of functional data exploring its most important features: such as depth measurements or functional outliers detection, among others.
It also helps to explain and model the relationship between a dependent variable and independent (regression models) and make predictions. Methods for supervised or unsupervised classification of a set of functional data regarding a feature of the data are also included. Finally, it can perform analysis of variance model (ANOVA) for functional data.

Details

Sections of fda.usc-package:

A.- Functional Data Representation
The functions included in this section allow to define, transform, manipulate and represent a functional dataset in many ways including derivatives, non-parametric kernel methods or basis representation.

B.- Functional Depth and Functional Outlier Detection

The functional data depth calculated by the different depth functions implemented that could be use as a measure of centrality or outlyingness.

B.1-Depth methods Depth:

B.2-Functional Outliers detection methods:

C.- Functional Regression Models

C.1. Functional explanatory covariate and scalar response
The functions included in this section allow the estimation of different functional regression models with a scalar response and a single functional explicative covariate.

C.2. Test for the functional linear model (FLM) with scalar response.

C.3. Functional and non functional explanatory covariates.
The functions in this section extends those regression models in previous section in several ways.

C.4. Functional response model (fregre.basis.fr) allows the estimation of functional regression models with a functional response and a single functional explicative covariate.

C.5. fregre.gls fits functional linear model using generalized least squares. fregre.igls fits iteratively a functional linear model using generalized least squares.

C.6. fregre.gsam.vs, Variable Selection using Functional Additive Models

D.- Functional Supervised Classification
This section allows the estimation of the groups in a training set of functional data fdata class by different nonparametric methods of supervised classification. Once these classifiers have been trained, they can be used to predict on new functional data.

Package allows the estimation of the groups in a training set of functional data by different methods of supervised classification.

D.1 Univariate predictor (x,y arguments, fdata class)

D.2 Multiple predictors (formula,data arguments, ldata class)

D.3 Depth classifiers (fdata or ldata class)

D.4 Functional Classification usign k-fold CV

E.- Functional Non-Supervised Classification
This section allows the estimation of the groups in a functional data set fdata class by kmeans method.

F.- Functional ANOVA

G.- Utilities and auxiliary functions:

Author(s)

Authors: Manuel Febrero Bande manuel.febrero@usc.es and Manuel Oviedo de la Fuente manuel.oviedo@udc.es

Contributors: Pedro Galeano, Alicia Nieto-Reyes, Eduardo Garcia-Portugues eduardo.garcia@usc.es and STAPH group https://www.math.univ-toulouse.fr/~ferraty/

Maintainer: Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Febrero-Bande, M., Oviedo de la Fuente, M. (2012). Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28. doi:10.18637/jss.v051.i04

Performance measures for regression and classification models

Description

cat2meas and tab2meas calculate the measures for a multiclass classification model.
pred2meas calculates the measures for a regression model.

Usage

cat2meas(yobs, ypred, measure = "accuracy", cost = rep(1, nlevels(yobs)))

tab2meas(tab, measure = "accuracy", cost = rep(1, nrow(tab)))

pred.MSE(yobs, ypred)

pred.RMSE(yobs, ypred)

pred.MAE(yobs, ypred)

pred2meas(yobs, ypred, measure = "RMSE")

Arguments

Details

• cat2meas compute tab=table(yobs,ypred) and calls tab2meas function.

• tab2meas function computes the following measures (see measure argument) for a binary classification model:

  • accuracy: Proportion of correct predictions. \frac{TP + TN}{TP + TN + FP + FN}

  • sensitivity, TPrate, recall: True Positive Rate or recall. \frac{TP}{TP + FN}

  • precision: Positive Predictive Value. \frac{TP}{TP + FP}

  • specificity, TNrate: True Negative Rate. \frac{TN}{TN + FP}

  • FPrate: False Positive Rate. \frac{FP}{TN + FP}

  • FNrate: False Negative Rate. \frac{FN}{TP + FN}

  • Fmeasure: Harmonic mean of precision and recall. \frac{2}{\frac{1}{\text{recall}} + \frac{1}{\text{precision}}}

  • Gmean: Geometric Mean of recall and specificity. \sqrt{\left(\frac{TP}{TP + FN}\right) \cdot \left(\frac{TN}{TN + FP}\right)}

  • kappa: Cohen's Kappa index. Kappa = \frac{P_o - P_e}{1 - P_e} where P_o is the proportion of observed agreement, \frac{TP + TN}{TP + TN + FP + FN}, and P_e is the proportion of agreement expected by chance, \frac{(TP + FP)(TP + FN) + (TN + FN)(TN + FP)}{(TP + TN + FP + FN)^2}.

  • cost: Weighted accuracy, calculated as \frac{\sum (\text{diag(tab)} / \text{rowSums(tab)} \cdot \text{cost})}{\sum(\text{cost})}

  • IOU: Mean Intersection over Union. \frac{TP}{TP + FN + FP}

  • IOU4class: Intersection over Union by class level. \frac{TP}{TP + FN + FP}#'

• pred2meas function computes the following measures of error, usign the measure argument, for observed and predicted vectors:

  • MSE: Mean squared error, \frac{\sum{(ypred- yobs)^2}}{n} .

  • RMSE: Root mean squared error \sqrt{\frac{\sum{(ypred- yobs)^2}}{n} }.

  • MAE: Mean Absolute Error, \frac{\sum |yobs - ypred|}{n}.

See Also

Other performance: weights4class()

aemet data

Description

Series of daily summaries of 73 spanish weather stations selected for the period 1980-2009. The dataset contains geographic information of each station and the average for the period 1980-2009 of daily temperature, daily precipitation and daily wind speed.

Format

Elements of aemet:
..$df: Data frame with information of each wheather station:

• ind: Indicated weather station.

• name: Station Name. 36 marked UTF-8 strings.

• province:Province (region) of Spain. 36 marked UTF-8 strings

• altitude: Altitude of the station (in meters).

• year.ini: Start year.

• year.end: End year.

• longitude: x geographic coordinate of the station (in decimal degrees).

• latitude: y geographic coordinate of the station (in decimal degrees).

The functional variables:

• ...$temp: mean curve of the average daily temperature for the period 1980-2009 (in degrees Celsius, marked with UTF-8 string). In leap years temperatures for February 28 and 29 were averaged.

• ...$wind.speed: mean curve of the average daily wind speed for the period 1980-2009 (in m/s).

• ...$logprec: mean curve of the log precipitation for the period 1980-2009 (in log mm). Negligible precipitation (less than 1 tenth of mm) is replaced by 0.05 and no precipitation (0.0 mm) is replaced by 0.01. Then the logarithm is applied.

Details

Meteorological State Agency of Spain (AEMET), https://www.aemet.es/es/portada. Government of Spain.
It marks 36 UTF-8 string of names of stations and 3 UTF-8 string names of provinces through the function iconv.

Author(s)

Manuel Febrero Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

Source

The data were obtained from the FTP of AEMET in 2009.

Examples

## Not run: 
data(aemet)
names(aemet)
names(aemet$df)
lat <- ifelse(aemet$df$latitude<31,"red","blue")
plot(aemet,col=lat)

## End(Not run)

DD-Classifier Based on DD-plot

Description

Fits Nonparametric Classification Procedure Based on DD–plot (depth-versus-depth plot) for G dimensions (G=g\times h, g levels and p data depth).

Usage

classif.DD(
  group,
  fdataobj,
  depth = "FM",
  classif = "glm",
  w,
  par.classif = list(),
  par.depth = list(),
  control = list(verbose = FALSE, draw = TRUE, col = NULL, alpha = 0.25)
)

Arguments

Details

Make the group classification of a training dataset using DD-classifier estimation in the following steps.

1. The function computes the selected depth measure of the points in fdataobj w.r.t. a subsample of each g level group and p data dimension (G=g \times p). The user can be specify the parameters for depth function in par.depth.

   (i) Type of depth function from functional data, see Depth:

   • "FM": Fraiman and Muniz depth.

   • "mode": h-modal depth.

   • "RT": Random Tukey depth.

   • "RP": Random projection depth.

   • "RPD": Double random projection depth.

   (ii) Type of depth function from multivariate functional data, see depth.mfdata:

   • "FMp": Fraiman and Muniz depth with common support. Suppose that all p-fdata objects have the same support (same rangevals); see depth.FMp.

   • "modep": h-modal depth using a p-dimensional metric; see depth.modep.

   • "RPp": Random projection depth using a p-variate depth with the projections; see depth.RPp.

   If the procedure requires to compute a distance such as in "knn" or "np" classifier or "mode" depth, the user must use a proper distance function: metric.lp for functional data and metric.dist for multivariate data.

   (iii) Type of depth function from multivariate data, see Depth.Multivariate:

   • "SD": Simplicial depth (for bivariate data).

   • "HS": Half-space depth.

   • "MhD": Mahalanobis depth.

   • "RD": Random projections depth.

   • "LD": Likelihood depth.

2. The function calculates the misclassification rate based on data depth computed in step (1) using the following classifiers.

   • "MaxD": Maximum depth.

   • "DD1": Search for the best separating polynomial of degree 1.

   • "DD2": Search for the best separating polynomial of degree 2.

   • "DD3": Search for the best separating polynomial of degree 3.

   • "glm": Logistic regression computed using Generalized Linear Models; see classif.glm.

   • "gam": Logistic regression computed using Generalized Additive Models; see classif.gsam.

   • "lda": Linear Discriminant Analysis computed using lda.

   • "qda": Quadratic Discriminant Analysis computed using qda.

   • "knn": k-Nearest Neighbour classification computed using classif.knn.

   • "np": Non-parametric Kernel classifier computed using classif.np.

   The user can be specify the parameters for classifier function in par.classif such as the smoothing parameter par.classif[["h"]], if classif="np" or the k-Nearest Neighbour par.classif[["knn"]], if classif="knn".

   In the case of polynomial classifier ("DD1", "DD2" and "DD3") uses the original procedure proposed by Li et al. (2012), by defalut rotating the DD-plot (to exchange abscise and ordinate) using in par.classif argument rotate=TRUE. Notice that the maximum depth classifier can be considered as a particular case of DD1, fixing the slope with a value of 1 (par.classif=list(pol=1)).

   The number of possible different polynomials depends on the sample size n and increases polynomially with order k. In the case of g groups, so the procedure applies some multiple-start optimization scheme to save time:

   • Generate all combinations of the elements of \( n \) taken \( k \) at a time: g \times combn(N,k) candidate solutions. When this number exceeds nmax=10000, a random sample of 10000 combinations is selected.

   • Smooth the empirical loss with the logistic function: 1/(1+e^{-tx}). The classification rule is constructed by optimizing the best noptim combinations in this random sample (by default, noptim=1 and tt=50/range(depth values)). Note that Li et al. found that the optimization results become stable for t \in [50, 200] when the depth is standardized with an upper bound of 1.

   The original procedure (Li et al. (2012)) not need to try many initial polynomials (nmax=1000) and that the procedure optimize the best (noptim=1), but we recommended to repeat the last step for different solutions, as for example nmax=250 and noptim=25. User can change the parameters pol, rotate, nmax, noptim and tt in the argument par.classif.

   The classif.DD procedure extends to multi-class problems by incorporating the method of majority voting in the case of polynomial classifier and the method One vs the Rest in the logistic case ("glm" and "gam").

Value

• group.est: Estimated vector groups by classified method selected.

• misclassification: Probability of misclassification.

• prob.classification: Probability of correct classification by group level.

• dep: Data frame with the depth of the curves for functional data (or points for multivariate data) in fdataobj w.r.t. each group level.

• depth: Character vector specifying the type of depth functions used.

• par.depth: List of parameters for depth function.

• classif: Type of classifier used.

• par.classif: List of parameters for classif procedure.

• w: Optional case weights.

• fit: Fitted object by classif method using the depth as covariate.

Author(s)

This version was created by Manuel Oviedo de la Fuente and Manuel Febrero Bande and includes the original version for polynomial classifier created by Jun Li, Juan A. Cuesta-Albertos and Regina Y. Liu.

References

Cuesta-Albertos, J.A., Febrero-Bande, M. and Oviedo de la Fuente, M. The DDG-classifier in the functional setting, (2017). Test, 26(1), 119-142. DOI: doi:10.1007/s11749-016-0502-6.

See Also

See Also as predict.classif.DD

Examples

## Not run: 
# DD-classif for functional data
data(tecator)
ab <- tecator$absorp.fdata
ab1 <- fdata.deriv(ab, nderiv = 1)
ab2 <- fdata.deriv(ab, nderiv = 2)
gfat <- factor(as.numeric(tecator$y$Fat>=15))

# DD-classif for p=1 functional  data set
out01 <- classif.DD(gfat,ab,depth="mode",classif="np")
out02 <- classif.DD(gfat,ab2,depth="mode",classif="np")
# DD-plot in gray scale
ctrl <- list(draw=T,col=gray(c(0,.5)),alpha=.2)
out02bis <- classif.DD(gfat,ab2,depth="mode",classif="np",control=ctrl)

# 2 depth functions (same curves) 
ldat <- mfdata("ab" = ab, "ab2" = ab2)
out03 <- classif.DD(gfat,list(ab2,ab2),depth=c("RP","mode"),classif="np")
# DD-classif for p=2 functional data set
# Weighted version 
out04 <- classif.DD(gfat, ldat, depth="mode",
                    classif="np", w=c(0.5,0.5))
# Model version
out05 <- classif.DD(gfat,ldat,depth="mode",classif="np")
# Integrated version (for multivariate functional data)
out06 <- classif.DD(gfat,ldat,depth="modep",classif="np")

# DD-classif for multivariate data
data(iris)
group <- iris[,5]
x <- iris[,1:4]
out07 <- classif.DD(group,x,depth="LD",classif="lda")
summary(out07)
out08 <- classif.DD(group, list(x,x), depth=c("MhD","LD"),
                    classif="lda")
summary(out08)

# DD-classif for functional data: g levels 
data(phoneme)
mlearn <- phoneme[["learn"]]
glearn <- as.numeric(phoneme[["classlearn"]])-1
out09 <- classif.DD(glearn,mlearn,depth="FM",classif="glm")
out10 <- classif.DD(glearn,list(mlearn,mlearn),depth=c("FM","RP"),classif="glm")
summary(out09)
summary(out10)

## End(Not run)

Classifier from Functional Data

Description

Classification of functional data using maximum depth.

Usage

classif.depth(
  group,
  fdataobj,
  newfdataobj,
  depth = "RP",
  par.depth = list(),
  CV = "none"
)

Arguments

Value

• group.est: Vector of classes of train sample data.

• group.pred: Vector of classes of test sample data.

• prob.classification: Probability of correct classification by group.

• max.prob: Highest probability of correct classification.

• fdataobj: fdata class object.

• group: Factor of length n.

Author(s)

Febrero-Bande, M. and Oviedo de la Fuente, M.

References

Cuevas, A., Febrero-Bande, M. and Fraiman, R. (2007). Robust estimation and classification for functional data via projection-based depth notions. Computational Statistics 22, 3, 481-496.

Examples

## Not run: 
data(phoneme)
mlearn<-phoneme[["learn"]]
mtest<-phoneme[["test"]]
glearn<-phoneme[["classlearn"]]
gtest<-phoneme[["classtest"]]

a1<-classif.depth(glearn,mlearn,depth="RP")
table(a1$group.est,glearn)
a2<-classif.depth(glearn,mlearn,depth="RP",CV=TRUE)
a3<-classif.depth(glearn,mlearn,depth="RP",CV=FALSE)
a4<-classif.depth(glearn,mlearn,mtest,"RP")
a5<-classif.depth(glearn,mlearn,mtest,"RP",CV=TRUE)     
table(a5$group.est,glearn)
a6<-classif.depth(glearn,mlearn,mtest,"RP",CV=FALSE)
table(a6$group.est,glearn)

## End(Not run)

Classification Fitting Functional Generalized Kernel Additive Models

Description

Computes functional classification using functional explanatory variables using backfitting algorithm.

Usage

classif.gkam(
  formula,
  data,
  weights = "equal",
  family = binomial(),
  par.metric = NULL,
  par.np = NULL,
  offset = NULL,
  prob = 0.5,
  type = "1vsall",
  control = NULL,
  ...
)

Arguments

Details

The first item in the data list is called "df" and is a data frame with the response, as glm.
Functional covariates of class fdata are introduced in the following items in the data list.

Value

Return gam object plus:

• formula: formula.

• data: List that containing the variables in the model.

• group: Factor of length n.

• group.est: Estimated vector groups.

• prob.classification: Probability of correct classification by group.

• prob.group: Matrix of predicted class probabilities. For each functional point shows the probability of each possible group membership.

• max.prob: Highest probability of correct classification.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Febrero-Bande M. and Gonzalez-Manteiga W. (2012). Generalized Additive Models for Functional Data. TEST. Springer-Velag. doi:10.1007/s11749-012-0308-0

McCullagh and Nelder (1989), Generalized Linear Models 2nd ed. Chapman and Hall.

Opsomer J.D. and Ruppert D.(1997). Fitting a bivariate additive model by local polynomial regression.Annals of Statistics, 25, 186-211.

See Also

See Also as: fregre.gkam.
Alternative method: classif.glm.

Examples

## Not run:  
## Time-consuming: selection of 2 levels 
data(phoneme)
mlearn<-phoneme[["learn"]][1:150]
glearn<-factor(phoneme[["classlearn"]][1:150])
dataf<-data.frame(glearn)
dat=list("df"=dataf,"x"=mlearn)
a1<-classif.gkam(glearn~x,data=dat)
summary(a1)
mtest<-phoneme[["test"]][1:150]
gtest<-factor(phoneme[["classtest"]][1:150])
newdat<-list("x"=mtest)
p1<-predict(a1,newdat)
table(gtest,p1)

## End(Not run)  

Classification Fitting Functional Generalized Linear Models

Description

Computes functional classification using functional (and non functional) explanatory variables by basis representation.

The first item in the data list is called "df" and is a data frame with the response and non functional explanatory variables, as glm.

Functional covariates of class fdata or fd are introduced in the following items in the data list.
basis.x is a list of basis for represent each functional covariate. The basis object can be created by the function: create.pc.basis, pca.fd create.pc.basis, create.fdata.basis o create.basis.
basis.b is a list of basis for represent each functional beta parameter. If basis.x is a list of functional principal components basis (see create.pc.basis or pca.fd) the argument basis.b is ignored.

Usage

classif.glm(
  formula,
  data,
  family = binomial(),
  weights = "equal",
  basis.x = NULL,
  basis.b = NULL,
  type = "1vsall",
  prob = 0.5,
  CV = FALSE,
  ...
)

Arguments

Value

Return glm object plus:

• formula: formula.

• data: List that containing the variables in the model.

• group: Factor of length n.

• group.est: Estimated vector groups.

• prob.classification: Probability of correct classification by group.

• prob.group: Matrix of predicted class probabilities. For each functional point shows the probability of each possible group membership.

• max.prob: Highest probability of correct classification.

Note

If the formula only contains a non functional explanatory variables (multivariate covariates), the function compute a standard glm procedure.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Ramsay, James O., and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York.

McCullagh and Nelder (1989), Generalized Linear Models 2nd ed. Chapman and Hall.

Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S, New York: Springer. Regression for R. R News 1(2):20-25

See Also

See Also as: fregre.glm.
classif.gsam and classif.gkam.

Examples

## Not run: 
data(phoneme)
ldat <- ldata("df" = data.frame(y = phoneme[["classlearn"]]),
             "x" = phoneme[["learn"]])
a1 <- classif.glm(y ~ x, data = ldat)
summary(a1)
newldat <- ldata("df" = data.frame(y = phoneme[["classtest"]]),
                "x" = phoneme[["test"]])
p1 <- predict(a1,newldat)
table(newldat$df$y,p1)
sum(p1==newldat$df$y)/250

## End(Not run)

Classification Fitting Functional Generalized Additive Models

Description

Computes functional classification using functional (and non functional) explanatory variables by basis representation.

Usage

classif.gsam(
  formula,
  data,
  family = binomial(),
  weights = "equal",
  basis.x = NULL,
  CV = FALSE,
  prob = 0.5,
  type = "1vsall",
  ...
)

Arguments

Details

The first item in the data list is called "df" and is a data frame with the response and non functional explanatory variables, as glm.

Functional covariates of class fdata or fd are introduced in the following items in the data list.
basis.x is a list of basis for represent each functional covariate. The basis object can be created by the function: create.pc.basis, pca.fd create.pc.basis, create.fdata.basis o create.basis.

Value

Return gam object plus:

• formula: formula.

• data: List that containing the variables in the model.

• group: Factor of length n.

• group.est: Estimated vector groups

• prob.classification: Probability of correct classification by group.

• prob.group: Matrix of predicted class probabilities. For each functional point shows the probability of each possible group membership.

• max.prob: Highest probability of correct classification.

• type: Type of classification scheme: 1 vs all or majority voting.

Note

If the formula only contains a non functional explanatory variables (multivariate covariates), the function compute a standard glm procedure.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Ramsay, James O., and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York.

McCullagh and Nelder (1989), Generalized Linear Models 2nd ed. Chapman and Hall.

Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S, New York: Springer. Regression for R. R News 1(2):20-25

See Also

See Also as: fregre.gsam.
Alternative method: classif.np, classif.glm and classif.gkam.

Examples

## Not run: 
data(phoneme)
ldat <- ldata("df" = data.frame(y = phoneme[["classlearn"]]),
              "x" = phoneme[["learn"]])
              classifKgroups <- fda.usc:::classifKgroups
a1 <- classif.gsam( y ~ s(x,k=3),data=ldat)
summary(a1)
newldat <- ldata("df" = data.frame(y = phoneme[["classtest"]]),
              "x" = phoneme[["test"]])
p1 <- predict(a1,newldat)
table(newldat$df$y,p1)
sum(p1==newldat$df$y)/250

## End(Not run)

Variable Selection in Functional Data Classification

Description

Computes classification by selecting the functional (and non functional) explanatory variables.

Usage

classif.gsam.vs(
  data = list(),
  y,
  x,
  family = binomial(),
  weights = "equal",
  basis.x = NULL,
  basis.b = NULL,
  type = "1vsall",
  prob = 0.5,
  alpha = 0.05,
  dcor.min = 0.01,
  smooth = TRUE,
  measure = "accuracy",
  xydist,
  ...
)

Arguments

Value

Return the final fitted model (same result of the classsification method) plus:

• dcor, matrix with the values of distance correlation for each pontential covariate (by column) and the residual of the model in each step (by row).

• i.predictor, vector with 1 if the variable is selected, 0 otherwise.

• ipredictor, vector with the name of selected variables (in order of selection)

Note

Adapted version from the original method in repression: fregre.gsam.vs.

Author(s)

Febrero-Bande, M. and Oviedo de la Fuente, M.

References

Febrero-Bande, M., Gonz\'alez-Manteiga, W. and Oviedo de la Fuente, M. Variable selection in functional additive regression models, (2018). Computational Statistics, 1-19. DOI: doi:10.1007/s00180-018-0844-5

See Also

See Also as: classif.gsam.

Examples

## Not run: 
data(tecator)
x <- tecator$absorp.fdata
x1 <- fdata.deriv(x)
x2 <- fdata.deriv(x,nderiv=2)
y <- factor(ifelse(tecator$y$Fat<12,0,1))
xcat0 <- cut(rnorm(length(y)),4)
xcat1 <- cut(tecator$y$Protein,4)
xcat2 <- cut(tecator$y$Water,4)
ind <- 1:129
dat    <- data.frame("Fat"=y, x1$data, xcat1, xcat2)
ldat <- ldata("df"=dat[ind,],"x"=x[ind,],"x1"=x1[ind,],"x2"=x2[ind,])
# 3 functionals (x,x1,x2), 3 factors (xcat0, xcat1, xcat2)
# and 100 scalars (impact poitns of x1) 

res.gam <- classif.gsam(Fat~s(x),data=ldat)
summary(res.gam)

# Time consuming
res.gam.vs <- classif.gsam.vs("Fat",data=ldat)
summary(res.gam.vs)
res.gam.vs$i.predictor
res.gam.vs$ipredictor

# Prediction 
newldat <- ldata("df"=dat[-ind,],"x"=x[-ind,],
                "x1"=x1[-ind,],"x2"=x2[-ind,])
pred.gam <- predict(res.gam,newldat)                
pred.gam.vs <- predict(res.gam.vs,newldat)
cat2meas(newldat$df$Fat, pred.gam)
cat2meas(newldat$df$Fat, pred.gam.vs)

## End(Not run)

Functional Classification usign k-fold CV

Description

Computes Functional Classification using k-fold cross-validation

Usage

classif.kfold(
  formula,
  data,
  classif = "classif.glm",
  par.classif,
  kfold = 10,
  param.kfold = NULL,
  measure = "accuracy",
  cost,
  models = FALSE,
  verbose = FALSE
)

Arguments

Details

Parameters for k-fold cross validation:

1. Number of basis elements:

   • Data-driven basis such as Functional Principal Componetns (PC). No implemented for PLS basis yet.

   • Fixed basis (bspline, fourier, etc.).

   Option used in some classifiers such as classif.glm, classif.gsam, classif.svm, etc.

2. Bandwidth parameter. Option used in non-parametric classificiation models such as classif.np and classif.gkam.

Value

Best fitted model computed by the k-fold CV using the method indicated in the classif argument and also returns:

1. param.min, value of parameter (or parameters) selected by k-fold CV.

2. params.error, k-fold CV error for each parameter combination.

3. pred.kfold, predicted response computed by k-fold CV.

4. model, if TRUE, list of models for each parameter combination.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

Examples

## Not run: 
data(tecator)    
cutpoint <- 18
tecator$y$class <- factor(ifelse(tecator$y$Fat<cutpoint,0,1))
table(tecator$y$class )
x <- tecator[[1]]
x2 <- fdata.deriv(tecator[[1]],2)
data  <-   list("df"=tecator$y,x=x,x2=x2)
formula  <-   formula(class~x+x2)

# ex: default excution of classifier (no k-fold CV)
classif="classif.glm"
out.default <- classif.kfold(formula, data, classif = classif)
out.default
out.default$param.min
out.default$params.error
summary(out.default)

# ex: Number of PC basis elements selected by 10-fold CV
# Logistic classifier
kfold = 10
param.kfold   <-   list("x"=list("pc"=c(1:8)),"x2"=list("pc"=c(1:8)))
out.kfold1   <-   classif.kfold(formula, data, classif = classif,
                            kfold = kfold,param.kfold = param.kfold)
out.kfold1$param.min
min(out.kfold1$params.error)
summary(out.kfold1)

# ex: Number of PC basis elements selected by 10-fold CV
# Logistic classifier with inverse weighting
out.kfold2 <- classif.kfold(formula, data, classif = classif,
                            par.classif=list("weights"="inverse"),
                            kfold = kfold,param.kfold = param.kfold)
out.kfold2$param.min
min(out.kfold2$params.error)
summary(out.kfold2)

# ex: Number of fourier  basis elements selected by 10-fold CV
# Logistic classifier 
ibase = seq(5,15,by=2)
param.kfold <- list("x"=list("fourier"=ibase),
                    "x2"=list("fourier"=ibase))
out.kfold3 <- classif.kfold(formula, data, classif = classif,
                            kfold = kfold,param.kfold = param.kfold)
out.kfold3$param.min
min(out.kfold3$params.error)
summary(out.kfold3)

# ex: Number of k-nearest neighbors selected by 10-fold CV
# non-parametric classifier  (only for a functional covariate)

output <- classif.kfold( class ~ x, data, classif = "classif.knn",
                       param.kfold= list("x"=list("knn"=c(3,5,9,13))))
output$param.min
output$params.error

output <- classif.kfold( class ~ x2, data, classif = "classif.knn",
                       param.kfold= list("x2"=list("knn"=c(3,5,9,13))))
output$param.min
output$params.error 

## End(Not run)
 

Functional classification using ML algotithms

Description

Computes functional classification using functional (and non functional) explanatory variables by rpart, nnet, svm or random forest model

Usage

classif.nnet(formula, data, basis.x = NULL, weights = "equal", size, ...)

classif.rpart(
  formula,
  data,
  basis.x = NULL,
  weights = "equal",
  type = "1vsall",
  ...
)

classif.svm(
  formula,
  data,
  basis.x = NULL,
  weights = "equal",
  type = "1vsall",
  ...
)

classif.ksvm(formula, data, basis.x = NULL, weights = "equal", ...)

classif.randomForest(
  formula,
  data,
  basis.x = NULL,
  weights = "equal",
  type = "1vsall",
  ...
)

classif.lda(
  formula,
  data,
  basis.x = NULL,
  weights = "equal",
  type = "1vsall",
  ...
)

classif.qda(
  formula,
  data,
  basis.x = NULL,
  weights = "equal",
  type = "1vsall",
  ...
)

classif.naiveBayes(formula, data, basis.x = NULL, laplace = 0, ...)

classif.cv.glmnet(formula, data, basis.x = NULL, weights = "equal", ...)

classif.gbm(formula, data, basis.x = NULL, weights = "equal", ...)

Arguments

Details

The first item in the data list is called "df" and is a data frame with the response and non functional explanatory variables, as glm.

Functional covariates of class fdata or fd are introduced in the following items in the data list.
basis.x is a list of basis for represent each functional covariate. The b object can be created by the function: create.pc.basis, pca.fd create.pc.basis, create.fdata.basis o create.basis.
basis.b is a list of basis for represent each functional beta parameter. If basis.x is a list of functional principal components basis (see create.pc.basis or pca.fd) the argument basis.b is ignored.

Value

Return classif object plus:

• formula: formula.

• data: List that containing the variables in the model.

• group: Factor of length n.

• group.est: Estimated vector groups.

• prob.classification: Probability of correct classification by group.

• prob.group: Matrix of predicted class probabilities. For each functional point shows the probability of each possible group membership.

• max.prob: Highest probability of correct classification.

• type: Type of classification scheme: 1 vs all or majority voting.

• fit: list of binary classification fitted models.

Note

Wrapper versions for multivariate and functional classification:

• classif.lda,classif.qda: uses lda and qda functions and requires MASS package.

• classif.nnet: uses nnet function and requires nnet package.

• classif.rpart: uses nnet function and requires rpart package.

• classif.svm, classif.naiveBayes: uses svm and naiveBayes functions and requires e1071 package.

• classif.ksvm: uses weighted.ksvm function and requires personalized package.

• classif.randomForest: uses randomForest function and requires randomForest package.

• classif.cv.glmnet: uses cv.glmnet function and requires glmnet package.

• classif.gbm: uses gbm function and requires gbm package.

Author(s)

Febrero-Bande, M. and Oviedo de la Fuente, M.

References

Ramsay, James O., and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York.

McCullagh and Nelder (1989), Generalized Linear Models 2nd ed. Chapman and Hall.

Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S, New York: Springer. Regression for R. R News 1(2):20-25

See Also

See Also as: rpart.
Alternative method: classif.np, classif.glm, classif.gsam and classif.gkam.

Examples

## Not run: 
data(phoneme)
mlearn<-phoneme[["learn"]]
glearn<-phoneme[["classlearn"]]
mtest<-phoneme[["test"]]
gtest<-phoneme[["classtest"]]
dataf<-data.frame(glearn)
dat=ldata("df"=dataf,"x"=mlearn)
a1<-classif.rpart(glearn~x,data=dat)
a2<-classif.nnet(glearn~x,data=dat)
a3<-classif.gbm(glearn~x,data=dat)
a4<-classif.randomForest(glearn~x,data=dat)
a5<-classif.cv.glmnet(glearn~x,data=dat)
newdat<-list("x"=mtest)
p1<-predict(a1,newdat,type="class")
p2<-predict(a2,newdat,type="class")
p3<-predict(a3,newdat,type="class")
p4<-predict(a4,newdat,type="class")
p5<-predict(a5,newdat,type="class")
mean(p1==gtest);mean(p2==gtest);mean(p3==gtest)
mean(p4==gtest);mean(p5==gtest)

## End(Not run)

Kernel Classifier from Functional Data

Description

Fits Nonparametric Supervised Classification for Functional Data.

Usage

classif.np(
  group,
  fdataobj,
  h = NULL,
  Ker = AKer.norm,
  metric,
  weights = "equal",
  type.S = S.NW,
  par.S = list(),
  ...
)

classif.knn(
  group,
  fdataobj,
  knn = NULL,
  metric,
  weights = "equal",
  par.S = list(),
  ...
)

classif.kernel(
  group,
  fdataobj,
  h = NULL,
  Ker = AKer.norm,
  metric,
  weights = "equal",
  par.S = list(),
  ...
)

Arguments

Details

Make the group classification of a training dataset using kernel or KNN estimation: Kernel.

Different types of metric funtions can be used.

Value

• fdataobj: fdata class object.

• group: Factor of length n.

• group.est: Estimated vector groups.

• prob.group: Matrix of predicted class probabilities. For each functional point shows the probability of each possible group membership.

• max.prob: Highest probability of correct classification.

• h.opt: Optimal smoothing parameter or bandwidht estimated.

• D: Matrix of distances of the optimal quantile distance hh.opt.

• prob.classification: Probability of correct classification by group.

• misclassification: Vector of probability of misclassification by number of neighbors knn.

• h: Vector of smoothing parameter or bandwidht.

• C: A call of function classif.kernel.

Note

If fdataobj is a data.frame the function considers the case of multivariate covariates.
metric.dist function is used to compute the distances between the rows of a data matrix (as dist function.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Ferraty, F. and Vieu, P. (2006). Nonparametric functional data analysis. Springer Series in Statistics, New York.

Ferraty, F. and Vieu, P. (2006). NPFDA in practice. Free access on line at https://www.math.univ-toulouse.fr/~ferraty/SOFTWARES/NPFDA/

See Also

See Also as predict.classif

Examples

## Not run: 
data(phoneme)
mlearn <- phoneme[["learn"]]
glearn <- phoneme[["classlearn"]]
h <- 9:19
out <- classif.np(glearn,mlearn,h=h)
summary(out)
head(round(out$prob.group,4))

## End(Not run)

Conditional Distribution Function

Description

Calculate the conditional distribution function of a scalar response with functional data.

Usage

cond.F(
  fdata0,
  y0,
  fdataobj,
  y,
  h = 0.15,
  g = 0.15,
  metric = metric.lp,
  Ker = list(AKer = AKer.epa, IKer = IKer.epa),
  ...
)

Arguments

Details

If x.dist=NULL the distance matrix between fdata objects is calculated by function passed in metric argument.

Value

• Fc: Conditional distribution function.

• y0: Vector of conditional response.

• g: Smoothing parameter or bandwidth of explanatory functional data (fdataobj).

• h: Smoothing parameter or bandwidth of respone, y.

• x.dist: Distance matrix between curves of fdataobj object.

• xy.dist: Distance matrix between cuves of fdataobj and fdata0 objects.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Ferraty, F. and Vieu, P. (2006). Nonparametric functional data analysis. Springer Series in Statistics, New York.

See Also

See Also as: cond.mode and cond.quantile.

Examples

## Not run: 
# Read data
n= 500
t= seq(0,1,len=101)
beta = t*sin(2*pi*t)^2
x = matrix(NA, ncol=101, nrow=n)
y=numeric(n)
x0<-rproc2fdata(n,seq(0,1,len=101),sigma="wiener")
x1<-rproc2fdata(n,seq(0,1,len=101),sigma=0.1)
x<-x0*3+x1
fbeta = fdata(beta,t)
y<-inprod.fdata(x,fbeta)+rnorm(n,sd=0.1) 
prx=x[1:100];pry=y[1:100]
ind=101;ind2=102:110
pr0=x[ind];pr10=x[ind2,]
ndist=61
gridy=seq(-1.598069,1.598069, len=ndist)

# Conditional Function
res1 = cond.F(pr10, gridy, prx, pry,p=1)
res2 = cond.F(pr10, gridy, prx, pry,h=0.3)
res3 = cond.F(pr10, gridy, prx, pry,g=0.25,h=0.3)

plot(res1$Fc[,1],type="l",ylim=c(0,1))
lines(res2$Fc[,1],type="l",col=2)
lines(res3$Fc[,1],type="l",col=3)

## End(Not run)

Conditional mode

Description

Computes the mode for conditional distribution function.

Usage

cond.mode(Fc, method = "monoH.FC", draw = TRUE)

Arguments

Details

The conditional mode is calculated as the maximum argument of the derivative of the conditional distribution function (density function f).

Value

Return the mode for conditional distribution function.

• mode.cond: Conditional mode.

• x: A grid of length n where the conditional density function is evaluated.

• f: The conditional density function evaluated at x.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Ferraty, F. and Vieu, P. (2006). Nonparametric functional data analysis. Springer Series in Statistics, New York.

See Also

See Also as: cond.F, cond.quantile and splinefun .

Examples

## Not run: 
n= 500
t= seq(0,1,len=101)
beta = t*sin(2*pi*t)^2
x = matrix(NA, ncol=101, nrow=n)
y=numeric(n)
x0<-rproc2fdata(n,seq(0,1,len=101),sigma="wiener")
x1<-rproc2fdata(n,seq(0,1,len=101),sigma=0.1)
x<-x0*3+x1
fbeta = fdata(beta,t)
y<-inprod.fdata(x,fbeta)+rnorm(n,sd=0.1)
prx=x[1:100];pry=y[1:100]
ind=101;ind2=101:110
pr0=x[ind];pr10=x[ind2]
ndist=161
gridy=seq(-1.598069,1.598069, len=ndist)
# Conditional Function
I=5
# Time consuming
res = cond.F(pr10[I], gridy, prx, pry, h=1)
mcond=cond.mode(res)
mcond2=cond.mode(res,method="diff")

## End(Not run)

Conditional quantile

Description

Computes the quantile for conditional distribution function.

Usage

cond.quantile(
  qua = 0.5,
  fdata0,
  fdataobj,
  y,
  fn,
  a = min(y),
  b = max(y),
  tol = 10^floor(log10(max(y) - min(y)) - 3),
  iter.max = 100,
  ...
)

Arguments

Value

Return the quantile for conditional distribution function.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Ferraty, F. and Vieu, P. (2006). Nonparametric functional data analysis. Springer Series in Statistics, New York.

See Also

See Also as: cond.F and cond.mode.

Examples

## Not run: 
n= 100
t= seq(0,1,len=101)
beta = t*sin(2*pi*t)^2
x = matrix(NA, ncol=101, nrow=n)
y=numeric(n)
x0<-rproc2fdata(n,seq(0,1,len=101),sigma="wiener")
x1<-rproc2fdata(n,seq(0,1,len=101),sigma=0.1)
x<-x0*3+x1
fbeta = fdata(beta,t)
y<-inprod.fdata(x,fbeta)+rnorm(n,sd=0.1)

prx=x[1:50];pry=y[1:50]
ind=50+1;ind2=51:60
pr0=x[ind];pr10=x[ind2]
ndist=161
gridy=seq(-1.598069,1.598069, len=ndist)
ind4=5
y0 = gridy[ind4]

# Conditional median
med=cond.quantile(qua=0.5,fdata0=pr0,fdataobj=prx,y=pry,fn=cond.F,h=1)

# Conditional CI 95% conditional
lo=cond.quantile(qua=0.025,fdata0=pr0,fdataobj=prx,y=pry,fn=cond.F,h=1)
up=cond.quantile(qua=0.975,fdata0=pr0,fdataobj=prx,y=pry,fn=cond.F,h=1)
print(c(lo,med,up))

## End(Not run)

Create Basis Set for Functional Data of fdata class

Description

Compute basis for functional data.

Usage

create.fdata.basis(
  fdataobj,
  l = 1:5,
  maxl = max(l),
  type.basis = "bspline",
  rangeval = fdataobj$rangeval,
  class.out = "fd"
)

create.pc.basis(
  fdataobj,
  l = 1:5,
  norm = TRUE,
  basis = NULL,
  lambda = 0,
  P = c(0, 0, 1),
  ...
)

create.pls.basis(
  fdataobj,
  y,
  l = 1:5,
  norm = TRUE,
  lambda = 0,
  P = c(0, 0, 1),
  ...
)

create.raw.fdata(fdataobj, l = 1:nrow(fdataobj))

Arguments

Value

• basis: basis object.

• x: if TRUE the value of the rotated data (the centred data multiplied by the basis matrix) is returned.

• mean: functional mean of fdataobj.

• df: degree of freedom.

• type: type of basis.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Ramsay, James O. and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York.

N. Kraemer, A.-L. Boulsteix, and G. Tutz (2008). Penalized Partial Least Squares with Applications to B-Spline Transformations and Functional Data. Chemometrics and Intelligent Laboratory Systems, 94, 60 - 69. doi:10.1016/j.chemolab.2008.06.009

See Also

See Also as create.basis and fdata2pc.

Examples

## Not run: 
data(tecator)
basis.pc<-create.pc.basis(tecator$absorp.fdata,c(1,4,5))
plot(basis.pc$basis,col=1)
summary(basis.pc)
basis.pls<-create.pls.basis(tecator$absorp.fdata,y=tecator$y[,1],c(1,4,5))
summary(basis.pls)
plot(basis.pls$basis,col=2)
summary(basis.pls)

basis.fd<-create.fdata.basis(tecator$absorp.fdata,c(1,4,5),
type.basis="fourier")
plot(basis.pc$basis)
basis.fdata<-create.fdata.basis(tecator$absorp.fdata,c(1,4,5),
type.basis="fourier",class.out="fdata")
plot(basis.fd,col=2,lty=1)
lines(basis.fdata,col=3,lty=1)

## End(Not run)

The cross-validation (CV) score

Description

Compute the leave-one-out cross-validation score.

Usage

CV.S(y, S, W = NULL, trim = 0, draw = FALSE, metric = metric.lp, ...)

Arguments

Details

A.-If trim=0:

CV(h)=\frac{1}{n} \sum_{i=1}^{n}{\Bigg(\frac{y_i-r_{i}(x_i)}{(1-S_{ii})}\Bigg)^{2}w(x_{i})}

S_{ii} is the ith diagonal element of the smoothing matrix S.

B.-If trim>0:

CV(h)=\frac{1}{l} \sum_{i=1}^{l}{\Bigg(\frac{y_i-r_{i}(x_i)}{(1-S_{ii})}\Bigg)^{2}w(x_{i})}

S_{ii} is the ith diagonal element of the smoothing matrix S and l the index of (1-trim) curves with less error.

Value

Returns CV score calculated for input parameters.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Wasserman, L. All of Nonparametric Statistics. Springer Texts in Statistics, 2006.

See Also

See Also as optim.np
Alternative method: GCV.S

Examples

## Not run: 
data(tecator)
x<-tecator$absorp.fdata
np<-ncol(x)
tt<-1:np
S1 <- S.NW(tt,3,Ker.epa)
S2 <- S.LLR(tt,3,Ker.epa)
S3 <- S.NW(tt,5,Ker.epa)
S4 <- S.LLR(tt,5,Ker.epa)
cv1 <- CV.S(x, S1)
cv2 <- CV.S(x, S2)
cv3 <- CV.S(x, S3)
cv4 <- CV.S(x, S4)
cv5 <- CV.S(x, S4,trim=0.1,draw=TRUE)
cv1;cv2;cv3;cv4;cv5
S6 <- S.KNN(tt,1,Ker.unif,cv=TRUE)
S7 <- S.KNN(tt,5,Ker.unif,cv=TRUE)
cv6 <- CV.S(x, S6)
cv7 <- CV.S(x, S7)
cv6;cv7

## End(Not run)
 

Distance Correlation Statistic and t-Test

Description

Distance correlation t-test of multivariate and functional independence (wrapper functions of energy package).

Usage

dcor.xy(x, y, test = TRUE, metric.x, metric.y, par.metric.x, par.metric.y, n)

dcor.dist(D1, D2)

bcdcor.dist(D1, D2, n)

dcor.test(D1, D2, n)

Arguments

Details

These wrapper functions extend the functions of the energy package for multivariate data to functional data. Distance correlation is a measure of dependence between random vectors introduced by Szekely, Rizzo, and Bakirov (2007). dcor.xy performs a nonparametric t-test of multivariate or functional independence in high dimension. The distribution of the test statistic is approximately Student t with n(n-3)/2-1 degrees of freedom and for n \geq 10 the statistic is approximately distributed as standard normal. Wrapper function of energy:::dcor.ttest. The t statistic is a transformation of a bias corrected version of distance correlation (see SR 2013 for details). Large values (upper tail) of the t statistic are significant.
dcor.test similar to dcor.xy but only for distance matrix. dcor.dist compute distance correlation statistic. Wrapper function of energy::dcor but only for distance matrix bcdcor.dist compute bias corrected distance correlation statistic. Wrapper function of energy:::bcdcor but only for distance matrix.

Value

dcor.test returns a list with class htest containing

• method: description of test.

• statistic: observed value of the test statistic.

• parameter: degrees of freedom.

• estimate: bias corrected distance correlation bcdcor(x,y).

• p.value: p-value of the t-test.

• data.name: description of data.

dcor.xy returns the previous list with class htest and

• D1: the distance matrix of x.

• D2: the distance matrix of y.

dcor.dist returns the distance correlation statistic.

bcdcor.dist returns the bias corrected distance correlation statistic.

Author(s)

Manuel Oviedo de la Fuente manuel.oviedo@udc.es and Manuel Febrero Bande

References

Szekely, G.J. and Rizzo, M.L. (2013). The distance correlation t-test of independence in high dimension. Journal of Multivariate Analysis, Volume 117, pp. 193-213.

Szekely, G.J., Rizzo, M.L., and Bakirov, N.K. (2007), Measuring and Testing Dependence by Correlation of Distances, Annals of Statistics, Vol. 35 No. 6, pp. 2769-2794.

See Also

metric.lp amd metric.dist.

Examples

## Not run:  
x<-rproc2fdata(100,1:50)
y<-rproc2fdata(100,1:50)
dcor.xy(x, y,test=TRUE)
dx <- metric.lp(x)
dy <- metric.lp(y)
dcor.test(dx, dy)
bcdcor.dist(dx, dy)
dcor.xy(x, y,test=FALSE)
dcor.dist(dx, dy)

## End(Not run)

Computation of depth measures for functional data

Description

Several depth measures can be computed for functional data for descriptive or classification purposes.

Usage

depth.mode(
  fdataobj,
  fdataori = fdataobj,
  trim = 0.25,
  metric = metric.lp,
  h = NULL,
  scale = FALSE,
  draw = FALSE,
  ...
)

depth.RP(
  fdataobj,
  fdataori = fdataobj,
  trim = 0.25,
  nproj = 50,
  proj = "vexponential",
  dfunc = "TD1",
  par.dfunc = list(),
  scale = FALSE,
  draw = FALSE,
  ...
)

depth.RPD(
  fdataobj,
  fdataori = fdataobj,
  nproj = 20,
  proj = 1,
  deriv = c(0, 1),
  trim = 0.25,
  dfunc2 = mdepth.LD,
  method = "fmm",
  draw = FALSE,
  ...
)

depth.RT(
  fdataobj,
  fdataori = fdataobj,
  trim = 0.25,
  nproj = 10,
  proj = 1,
  xeps = 1e-07,
  draw = FALSE,
  ...
)

depth.KFSD(
  fdataobj,
  fdataori = fdataobj,
  trim = 0.25,
  h = NULL,
  scale = FALSE,
  draw = FALSE
)

depth.FSD(
  fdataobj,
  fdataori = fdataobj,
  trim = 0.25,
  scale = FALSE,
  draw = FALSE
)

depth.FM(
  fdataobj,
  fdataori = fdataobj,
  trim = 0.25,
  scale = FALSE,
  dfunc = "FM1",
  par.dfunc = list(scale = TRUE),
  draw = FALSE
)

Arguments

Details

Type of depth functions: Fraiman and Muniz (FM) depth, modal depth, random Tukey (RT), random projection (RP) depth and double random projection depth (RPD).

• depth.FM: computes the integration of an univariate depth along the axis x (see Fraiman and Muniz 2001). It is also known as Integrated Depth.

• depth.mode: implements the modal depth (see Cuevas et al 2007).

• depth.RT: implements the Random Tukey depth (see Cuesta–Albertos and Nieto–Reyes 2008).

• depth.RP: computes the Random Projection depth (see Cuevas et al. 2007).

• depth.RPD: implements a depth measure based on random projections possibly using several derivatives (see Cuevas et al. 2007).

• depth.FSD: computes the Functional Spatial Depth (see Sguera et al. 2014).

• depth.KFSD: implements the Kernelized Functional Spatial Depth (see Sguera et al. 2014).

• depth.mode: calculates the depth of a datum accounting the number of curves in its neighbourhood. By default, the distance is calculated using metric.lp function although any other distance could be employed through argument metric (with the general pattern USER.DIST(fdataobj,fdataori)).

• depth.RP: summarizes the random projections through averages whereas the depth.RT function uses the minimum of all projections.

• depth.RPD: involves the original trajectories and the derivatives of each curve in two steps. It builds random projections for the function and their derivatives (indicated in the parameter deriv) and then applies a depth function (by default depth.mode) to this set of random projections (by default the Tukey one).

• depth.FSD and depth.KFSD: are the implementations of the default versions of the functional spatial depths proposed in Sguera et al 2014. At this moment, it is not possible to change the kernel in the second one.

Value

Return a list with:

• median: Deepest curve.

• lmed: Index deepest element median.

• mtrim: fdata class object with the average from the (1-trim)% deepest curves.

• ltrim: Indexes of curves that conform the trimmed mean mtrim.

• dep: Depth of each curve of fdataobj w.r.t. fdataori.

• dep.ori: Depth of each curve of fdataori w.r.t. fdataori.

• proj: The projection value of each point on the curves.

• dist: Distance matrix between curves or functional data.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Cuevas, A., Febrero-Bande, M., Fraiman, R. (2007). Robust estimation and classification for functional data via projection-based depth notions. Computational Statistics 22, 3, 481-496.

Fraiman R, Muniz G. (2001). Trimmed means for functional data. Test 10: 419-440.

Cuesta–Albertos, JA, Nieto–Reyes, A. (2008) The Random Tukey Depth. Computational Statistics and Data Analysis Vol. 52, Issue 11, 4979-4988.

Febrero-Bande, M, Oviedo de la Fuente, M. (2012). Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28. https://www.jstatsoft.org/v51/i04/

Sguera C, Galeano P, Lillo R (2014). Spatial depth based classification for functional data. TEST 23(4):725–750.

See Also

See Also as Descriptive.

Examples

## Not run: 
#Ex: CanadianWeather data
tt=1:365
fdataobj<-fdata(t(CanadianWeather$dailyAv[,,1]),tt)
# Fraiman-Muniz Depth
out.FM=depth.FM(fdataobj,trim=0.1,draw=TRUE)
#Modal Depth
out.mode=depth.mode(fdataobj,trim=0.1,draw=TRUE)
out.RP=depth.RP(fdataobj,trim=0.1,draw=TRUE)
out.RT=depth.RT(fdataobj,trim=0.1,draw=TRUE)
out.FSD=depth.FSD(fdataobj,trim=0.1,draw=TRUE)
out.KFSD=depth.KFSD(fdataobj,trim=0.1,draw=TRUE)
## Double Random Projections
out.RPD=depth.RPD(fdataobj,deriv=c(0,1),dfunc2=mdepth.LD,
trim=0.1,draw=TRUE)
out<-c(out.FM$mtrim,out.mode$mtrim,out.RP$mtrim,out.RPD$mtrim)
plot(fdataobj,col="grey")
lines(out)
cdep<-cbind(out.FM$dep,out.mode$dep,out.RP$dep,out.RT$dep,out.FSD$dep,out.KFSD$dep)
colnames(cdep)<-c("FM","mode","RP","RT","FSD","KFSD")
pairs(cdep)
round(cor(cdep),2)

## End(Not run)

Provides the depth measure for multivariate data

Description

Compute measure of centrality of the multivariate data. Type of depth function: simplicial depth (SD), Mahalanobis depth (MhD), Random Half–Space depth (HS), random projection depth (RP) and Likelihood Depth (LD).

Usage

mdepth.LD(x, xx = x, metric = metric.dist, h = NULL, scale = FALSE, ...)

mdepth.HS(x, xx = x, proj = 50, scale = FALSE, xeps = 1e-15, random = FALSE)

mdepth.RP(x, xx = x, proj = 50, scale = FALSE)

mdepth.MhD(x, xx = x, scale = FALSE)

mdepth.KFSD(x, xx = x, trim = 0.25, h = NULL, scale = FALSE, draw = FALSE)

mdepth.FSD(x, xx = x, trim = 0.25, scale = FALSE, draw = FALSE)

mdepth.FM(x, xx = x, scale = FALSE, dfunc = "TD1")

mdepth.TD(x, xx = x, xeps = 1e-15, scale = FALSE)

mdepth.SD(x, xx = NULL, scale = FALSE)

Arguments

Details

Type of depth measures:

• The mdepth.SD calculates the simplicial depth (HD) of the points in x w.r.t. xx (for bivariate data).

• The mdepth.HS function calculates the random half–space depth (HS) of the points in x w.r.t. xx based on random projections proj.

• The mdepth.MhD function calculates the Mahalanobis depth (MhD) of the points in x w.r.t. xx.

• The mdepth.RP calculates the random' projection depth (RP) of the points in x w.r.t. xx based on random projections proj.

• The mdepth.LD calculates the Likelihood depth (LD) of the points in x w.r.t. xx.

• The mdepth.TD function provides the Tukey depth measure for multivariate data.

Value

• lmed: Index of the deepest element (median) of xx.

• ltrim: Index of the set of points x with trimmed mean mtrim.

• dep: Depth of each point in x with respect to xx.

• proj: The projection value of each point onto the set of points.

• x: A set of points to be evaluated.

• xx: A reference sample.

• name: Name of the depth method.

Author(s)

mdepth.RP, mdepth.MhD and mdepth.HS are versions created by Manuel Febrero Bande and Manuel Oviedo de la Fuente of the original version created by Jun Li, Juan A. Cuesta Albertos and Regina Y. Liu for polynomial classifier.

References

Liu, R. Y., Parelius, J. M., and Singh, K. (1999). Multivariate analysis by data depth: descriptive statistics, graphics and inference,(with discussion and a rejoinder by Liu and Singh). The Annals of Statistics, 27(3), 783-858.

See Also

Functional depth functions: depth.FM, depth.mode, depth.RP, depth.RPD and depth.RT.

Examples

## Not run: 
data(iris)
group<-iris[,5]
x<-iris[,1:2]
                                  
MhD<-mdepth.MhD(x)
PD<-mdepth.RP(x)
HD<-mdepth.HS(x)
SD<-mdepth.SD(x)

x.setosa<-x[group=="setosa",]
x.versicolor<-x[group=="versicolor",] 
x.virginica<-x[group=="virginica",]
d1<-mdepth.SD(x,x.setosa)$dep
d2<-mdepth.SD(x,x.versicolor)$dep
d3<-mdepth.SD(x,x.virginica)$dep

## End(Not run)

Provides the depth measure for a list of p–functional data objects

Description

This function computes the depth measure for a list of p–functional data objects. The procedure extends the Fraiman and Muniz (FM), modal, and random project depth functions from 1 functional dataset to p functional datasets.

Usage

depth.modep(
  mfdata,
  mfdataref = mfdata,
  h = NULL,
  metric,
  par.metric = list(),
  method = "euclidean",
  scale = FALSE,
  trim = 0.25,
  draw = FALSE,
  ask = FALSE
)

depth.RPp(
  mfdata,
  mfdataref = mfdata,
  nproj = 50,
  proj = "vexponential",
  trim = 0.25,
  dfunc = "mdepth.TD",
  par.dfunc = list(scale = TRUE),
  draw = FALSE,
  ask = FALSE
)

depth.FMp(
  mfdata,
  mfdataref = mfdata,
  trim = 0.25,
  dfunc = "mdepth.MhD",
  par.dfunc = list(scale = FALSE),
  draw = FALSE,
  ask = FALSE,
  ...
)

Arguments

Details

• depth.FMp, this procedure suposes that each curve of the mfdataobj have the same support [0,T] (same argvals and rangeval). The FMp depth is defined as: FM_i^p =\int_{0}^{T}Z_i^p(t)dt where Z_i^p(t) is a p–variate depth of the vector (x_i^1(t),\ldots,x_i^p(t)) w.r.t. the sample at t. derivatives. In this case,note solo un dato funcional se reduce depth.FM=depth.FM1

• The depth.RPp function calculates the depth in two steps. It builds random projections for the each curve of the mfdata w.r.t. each curve of the mfdataref object. Then it applyes a multivariate depth function specified in dfunc argument to the set of random projections. This procedure is a generalization of Random Projection with derivatives (RPD) implemented in depth.RPD function. Now, the procedure computes a p-variate depth with the projections using the p functional dataset.

• The modal depth depth.modep function calculates the depth in three steps. First, the function calculates a suitable metrics or semi–metrics m_1+\cdots+m_p for each curve of the mfdata w.r.t. each curve in the mfdataref object using the metric and par.metric arguments, see metric.lp or semimetric.NPFDA for more details. Second, the function uses the p–dimensional metrics to construct a new metric, specified in method argument, by default if method="euclidean", i.e. m:=\sqrt{m_1^2+\cdots+m_p^2}. Finally, the empirical h–depth is computed as:

  \hat{f}_h(x_0)=N^{-1}\sum_{i=1}^{N}{K(m/h)}

  where x is dataset with p observed fucntional data, m is a suitable metric or semi–metric, K(t) is an asymmetric kernel function and h is the bandwidth parameter.

Value

• lmed: Index deepest element median.

• ltrim: Index of curves with trimmed mean mtrim.

• dep: Depth of each curve of fdataobj w.r.t. fdataori.

• dfunc: Second depth function used as multivariate depth, see details section.

• par.dfunc: List of parameters for the dfunc depth function.

• proj: The projection value of each point on the curves.

• dist: Distance matrix between curves or functional data.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Cuevas, A., Febrero-Bande, M. and Fraiman, R. (2007). Robust estimation and classification for functional data via projection-based depth notions. Computational Statistics 22, 3, 481-496.

See Also

See Also as Descriptive.

Examples

## Not run: 
data(tecator)
xx<-tecator$absorp
xx1<-fdata.deriv(xx,1)
lx<-list(xx=xx,xx=xx1)
# Fraiman-Muniz Depth
par.df<-list(scale =TRUE)
out.FM1p=depth.FMp(lx,trim=0.1,draw=TRUE, par.dfunc = par.df)
out.FM2p=depth.FMp(lx,trim=0.1,dfunc="mdepth.LD",
par.dfunc = par.df, draw=TRUE)

# Random Project Depth
out.RP1p=depth.RPp(lx,trim=0.1,dfunc="mdepth.TD",
draw=TRUE,par.dfunc = par.df)
out.RP2p=depth.RPp(lx,trim=0.1,dfunc="mdepth.LD",
draw=TRUE,par.dfunc = par.df)

#Modal Depth
out.mode1p=depth.modep(lx,trim=0.1,draw=T,scale=T)
out.mode2p=depth.modep(lx,trim=0.1,method="manhattan",
draw=T,scale=T)

par(mfrow=c(2,3))
plot(out.FM1p$dep,out.FM2p$dep)
plot(out.RP1p$dep,out.RP2p$dep)
plot(out.mode1p$dep,out.mode2p$dep)
plot(out.FM1p$dep,out.RP1p$dep)
plot(out.RP1p$dep,out.mode1p$dep)
plot(out.FM1p$dep,out.mode1p$dep)

## End(Not run)

Descriptive measures for functional data.

Description

Central and dispersion measures for functional data.

Usage

func.mean(x)

func.var(fdataobj)

func.trim.FM(fdataobj, ...)

func.trim.mode(fdataobj, ...)

func.trim.RP(fdataobj, ...)

func.trim.RT(fdataobj, ...)

func.trim.RPD(fdataobj, ...)

func.med.FM(fdataobj, ...)

func.med.mode(fdataobj, ...)

func.med.RP(fdataobj, ...)

func.med.RT(fdataobj, ...)

func.med.RPD(fdataobj, ...)

func.trimvar.FM(fdataobj, ...)

func.trimvar.mode(fdataobj, ...)

func.trimvar.RP(fdataobj, ...)

func.trimvar.RPD(fdataobj, ...)

func.trim.RT(fdataobj, ...)

func.med.RT(fdataobj, ...)

func.trimvar.RT(fdataobj, ...)

func.mean.formula(formula, data = NULL, ..., drop = FALSE)

Arguments

Value

func.mean.formula The value returned from split is a list of fdata containing the mean curves
for the groups. The components of the list are named by the levels of f (after converting to a factor, or if already a factor and drop = TRUE, dropping unused levels).

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Febrero-Bande, M., Oviedo de la Fuente, M. (2012). Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28. https://www.jstatsoft.org/v51/i04/

Examples

## Not run: 
#' # Example with Montreal Daily Temperature (fda-package)
fdataobj<-fdata(MontrealTemp)

# Measures of central tendency by group
fac<-factor(c(rep(1,len=17),rep(2,len=17)))
ldata=list("df"=data.frame(fac),"fdataobj"=fdataobj)
a1<-func.mean.formula(fdataobj~fac,data=ldata)
plot(a1)

# Measures of central tendency
a1<-func.mean(fdataobj)
a2<-func.trim.FM(fdataobj)
a3<-func.trim.mode(fdataobj)
a4<-func.trim.RP(fdataobj)
# a5<-func.trim.RPD(fdataobj,deriv=c(0,1)) # Time-consuming
a6<-func.med.FM(fdataobj)
a7<-func.med.mode(fdataobj)
a8<-func.med.RP(fdataobj)
# a9<-func.med.RPD(fdataobj,deriv=c(0,1)) # Time-consuming
# a10<-func.med.RT(fdataobj)

par(mfrow=c(1,2))
plot(c(a1,a2,a3,a4),ylim=c(-26,29),main="Central tendency: trimmed mean")
plot(c(a1,a6,a7,a8),ylim=c(-26,29),main="Central tendency: median")

## Measures of dispersion
b1<-func.var(fdataobj)
b2<-func.trimvar.FM(fdataobj)
b3<-func.trimvar.FM(fdataobj,trim=0.1)
b4<-func.trimvar.mode(fdataobj)
b5<-func.trimvar.mode(fdataobj,p=1)
b6<-func.trimvar.RP(fdataobj)
b7<-func.trimvar.RPD(fdataobj)
b8<-func.trimvar.RPD(fdataobj)
b9<-func.trimvar.RPD(fdataobj,deriv=c(0,1))
dev.new()
par(mfrow=c(1,2))
plot(c(b1,b2,b3,b4,b5),ylim=c(0,79),main="Measures of dispersion I")
plot(c(b1,b6,b7,b8,b9),ylim=c(0,79),main="Measures of dispersion II")

## End(Not run)

The deviance score

Description

Returns the deviance of a fitted model object by GCV score.

Usage

dev.S(
  y,
  S,
  obs,
  family = gaussian(),
  off,
  offdf,
  criteria = "GCV",
  W = diag(1, ncol = ncol(S), nrow = nrow(S)),
  trim = 0,
  draw = FALSE,
  ...
)

Arguments

Details

Up to a constant, minus twice the maximized log-likelihood. Where sensible, the constant is chosen so that a saturated model has deviance zero.

GCV(h)=p(h) \Xi(n^{-1}h^{-1})

Where

p(h)=\frac{1}{n} \sum_{i=1}^{n}{\Big(y_i-r_{i}(x_i)\Big)^{2}w(x_i)}

and penalty function

\Xi()

can be selected from the following criteria:

Generalized Cross-validation (GCV):

\Xi_{GCV}(n^{-1}h^{-1})=(1-n^{-1}S_{ii})^{-2}

Akaike's Information Criterion (AIC):

\Xi_{AIC}(n^{-1}h^{-1})=exp(2n^{-1}S_{ii})

Finite Prediction Error (FPE)

\Xi_{FPE}(n^{-1}h^{-1})=\frac{(1+n^{-1}S_{ii})}{(1-n^{-1}S_{ii})}

Shibata's model selector (Shibata):

\Xi_{Shibata}(n^{-1}h^{-1})=(1+2n^{-1}S_{ii})

Rice's bandwidth selector (Rice):

\Xi_{Rice}(n^{-1}h^{-1})=(1-2n^{-1}S_{ii})^{-1}

Value

Returns GCV score calculated for input parameters.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Wasserman, L. All of Nonparametric Statistics. Springer Texts in Statistics, 2006.

Hardle, W. Applied Nonparametric Regression. Cambridge University Press, 1994.

Febrero-Bande, M., Oviedo de la Fuente, M. (2012). Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28. https://www.jstatsoft.org/v51/i04/

See Also

See Also as GCV.S.
Alternative method: CV.S

Examples

data(phoneme)
mlearn<-phoneme$learn
np<-ncol(mlearn)
tt<-mlearn[["argvals"]]
S1 <- S.NW(tt,2.5)
gcv1 <- dev.S(mlearn$data[1,],obs=(sample(150)), 
S1,off=rep(1,150),offdf=3)
gcv2 <- dev.S(mlearn$data[1,],obs=sort(sample(150)), 
S1,off=rep(1,150),offdf=3)

Delsol, Ferraty and Vieu test for no functional-scalar interaction

Description

The function dfv.test tests the null hypothesis of no interaction between a functional covariate and a scalar response in a general framework. The null hypothesis is

H_0:\,m(X)=0,

where m(\cdot) denotes the regression function of the functional variate X over the centred scalar response Y (E[Y]=0). The null hypothesis is tested by the smoothed integrated square error of the response (see Details).

Usage

dfv.statistic(
  X.fdata,
  Y,
  h = quantile(x = metric.lp(X.fdata), probs = c(0.05, 0.1, 0.15, 0.25, 0.5)),
  K = function(x) 2 * dnorm(abs(x)),
  weights = rep(1, dim(X.fdata$data)[1]),
  d = metric.lp,
  dist = NULL
)

dfv.test(
  X.fdata,
  Y,
  B = 5000,
  h = quantile(x = metric.lp(X.fdata), probs = c(0.05, 0.1, 0.15, 0.25, 0.5)),
  K = function(x) 2 * dnorm(abs(x)),
  weights = rep(1, dim(X.fdata$data)[1]),
  d = metric.lp,
  verbose = TRUE
)

Arguments

Details

The Delsol, Ferraty and Vieu statistic is defined as

T_n=\int\bigg(\sum_{i=1}^n(Y_i-m(X_i))K\bigg(\frac{d(X,X_i)}{h}\bigg)\bigg)^2\omega(X)dP_X(X)

and in the case of no interaction with centred scalar response (when H_0:\,m(X)=0 holds), its sample version is computed from

T_n=\frac{1}{n}\sum_{j=1}^n\bigg(\sum_{i=1}^n Y_iK\bigg(\frac{d(X_j,X_i)}{h}\bigg)\bigg)^2\omega(X_j).

The sample version implemented here does not consider a splitting of the sample, as the authors comment in their paper. The statistic is computed by the function dfv.statistic and, before applying the test, the response Y is centred. The distribution of the test statistic is approximated by a wild bootstrap on the residuals, using the golden section bootstrap.

Please note that if a grid of bandwidths is passed, a harmless warning message will prompt at the end of the test (it comes from returning several p-values in the htest class).

Value

The value of dfv.statistic is a vector of length length(h) with the values of the statistic for each bandwidth. The value of dfv.test is an object with class "htest" whose underlying structure is a list containing the following components:

• statistic: The value of the Delsol, Ferraty and Vieu test statistic.

• boot.statistics: A vector of length B with the values of the bootstrap test statistics.

• p.value: The p-value of the test.

• method: The character string "Delsol, Ferraty and Vieu test for no functional-scalar interaction".

• B: The number of bootstrap replicates used.

• h: Bandwidth parameters for the test.

• K: Kernel function used.

• weights: The weights considered.

• d: Matrix of distances of the functional data.

• data.name: The character string "Y=0+e".

Note

No NA's are allowed neither in the functional covariate nor in the scalar response.

Author(s)

Eduardo Garcia-Portugues. Please, report bugs and suggestions to eduardo.garcia.portugues@uc3m.es

References

Delsol, L., Ferraty, F. and Vieu, P. (2011). Structural test in regression on functional variables. Journal of Multivariate Analysis, 102, 422-447. doi:10.1016/j.jmva.2010.10.003

Delsol, L. (2013). No effect tests in regression on functional variable and some applications to spectrometric studies. Computational Statistics, 28(4), 1775-1811. doi:10.1007/s00180-012-0378-1

See Also

rwild, flm.test, flm.Ftest, fregre.np

Examples

## Not run: 
## Simulated example ##
X=rproc2fdata(n=50,t=seq(0,1,l=101),sigma="OU")

beta0=fdata(mdata=rep(0,length=101)+rnorm(101,sd=0.05),
argvals=seq(0,1,l=101),rangeval=c(0,1))
beta1=fdata(mdata=cos(2*pi*seq(0,1,l=101))-(seq(0,1,l=101)-0.5)^2+
rnorm(101,sd=0.05),argvals=seq(0,1,l=101),rangeval=c(0,1))

# Null hypothesis holds
Y0=drop(inprod.fdata(X,beta0)+rnorm(50,sd=0.1))

# Null hypothesis does not hold
Y1=drop(inprod.fdata(X,beta1)+rnorm(50,sd=0.1))

# We use the CV bandwidth given by fregre.np
# Do not reject H0
dfv.test(X,Y0,h=fregre.np(X,Y0)$h.opt,B=100)
# dfv.test(X,Y0,B=5000)

# Reject H0
dfv.test(X,Y1,B=100)
# dfv.test(X,Y1,B=5000)

## End(Not run)

Proximities between functional data

Description

Computes the cosine correlation distance between two functional dataset of class fdata.

Usage

dis.cos.cor(fdata1, fdata2 = NULL, as.dis = FALSE)

Arguments

Value

Returns a proximity/distance matrix (depending on as.dis) between functional data.

References

Kemmeren P, van Berkum NL, Vilo J, et al. (2002). Protein Interaction Verification and Functional Annotation by Integrated Analysis of Genome-Scale Data . Mol Cell. 2002 9(5):1133-43.

See Also

See also metric.lp and semimetric.NPFDA

Examples

## Not run: 
 r1<-rnorm(1001,sd=.01)
 r2<-rnorm(1001,sd=.01)
 x<-seq(0,2*pi,length=1001)
 fx<-fdata(sin(x)/sqrt(pi)+r1,x)
 dis.cos.cor(fx,fx)
 dis.cos.cor(c(fx,fx),as.dis=TRUE)
 fx0<-fdata(rep(0,length(x))+r2,x)
 plot(c(fx,fx0))
 dis.cos.cor(c(fx,fx0),as.dis=TRUE)
 
## End(Not run)
 

ANOVA for heteroscedastic data

Description

Univariate ANOVA for heteroscedastic data.

Usage

fanova.hetero(object = NULL, formula, pr = FALSE, contrast = NULL, ...)

Arguments

Details

This function fits a univariate analysis of variance model and allows calculate special contrasts defined by the user. The list of special contrast to be used for some of the factors in the formula. Each matrix of the list has r rows and r-1 columns.

The user can also request special predetermined contrasts, for example using contr.helmert, contr.sum or contr.treatment functions.

Value

Return:

• ans: A list with components including: the Beta estimation Est, the factor degrees of freedom df1, the residual degrees of freedom df2, and the p-value for each factor.

• contrast: List of special contrasts.

Note

anova.hetero deprecated

It only works with categorical variables.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Brunner, E., Dette, H., Munk, A. Box-Type Approximations in Nonparametric Factorial Designs. Journal of the American Statistical Association, Vol. 92, No. 440 (Dec., 1997), pp. 1494-1502.

See Also

See Also as: fanova.RPm

Examples

## Not run: 
data(phoneme)
ind=1 # beetwen 1:150
fdataobj=data.frame(phoneme$learn[["data"]][,ind])
n=dim(fdataobj)[1]
group<-factor(phoneme$classlearn)

#ex 1: real factor and random factor
group.rand=as.factor(sample(rep(1:3,n),n))
f=data.frame(group,group.rand)
mm=data.frame(fdataobj,f)
colnames(mm)=c("value","group","group.rand")
out1=fanova.hetero(object=mm[,-2],value~group.rand,pr=FALSE)
out2=fanova.hetero(object=mm[,-3],value~group,pr=FALSE)
out1
out2

#ex 2: real factor, random factor and  special contrasts
cr5=contr.sum(5)  #each level vs last level
cr3=c(1,0,-1)			#first level vs last level
out.contrast=fanova.hetero(object=mm[,-3],value~group,pr=FALSE,
contrast=list(group=cr5))
out.contrast

## End(Not run)     

One–way anova model for functional data

Description

One–way anova model for k independent samples of functional data. The function contrasts the null hypothesis of equality of mean functions of functional data based on the an asymptotic version of the anova F–test.

H_0:\, m_1=\ldots=m_k

Usage

fanova.onefactor(
  object,
  group,
  nboot = 100,
  plot = FALSE,
  verbose = FALSE,
  ...
)

Arguments

Details

The function returns the p–value of test using one–way anova model over nboot runs.

Value

Returns:

• p-value: Probability of rejecting the null hypothesis H0 at a significance level.

• stat: Statistic value of the test.

• wm: Statistic values of bootstrap resamples.

Note

anova.onefactor deprecated.

Author(s)

Juan A. Cuesta-Albertos, Manuel Febrero-Bande, Manuel Oviedo de la Fuente
manuel.oviedo@udc.es

References

Cuevas, A., Febrero, M., & Fraiman, R. (2004). An anova test for functional data. Computational statistics & data analysis, 47(1), 111-122.

See Also

See Also as: fanova.RPm

Examples

## Not run: 
data(MCO)
grupo<-MCO$classintact
datos<-MCO$intact
res=fanova.onefactor(datos,grupo,nboot=50,plot=TRUE)
grupo <- MCO$classpermea
datos <- MCO$permea
res=fanova.onefactor(datos,grupo,nboot=50,plot=TRUE)

## End(Not run)  

Functional ANOVA with Random Project.

Description

The procedure is based on the analysis of randomly chosen one-dimensional projections. The function tests ANOVA models for functional data with continuous covariates and perform special contrasts for the factors in the formula.

Usage

fanova.RPm(
  object,
  formula,
  data.fac,
  RP = min(30, ncol(object)),
  alpha = 0.95,
  zproj = NULL,
  par.zproj = list(norm = TRUE),
  hetero = TRUE,
  pr = FALSE,
  w = rep(1, ncol(object)),
  nboot = 0,
  contrast = NULL,
  ...
)

## S3 method for class 'fanova.RPm'
summary(object, ndec = NULL, ...)

Arguments

Details

zproj allows to change the generator process of the projections. This can be done through the inclusion of a function or a collection of projections generated outside the function. By default, for a functional problem, the function rproc2fdata is used. For multivariate problems, if no function is included, the projections are generated by a normalized gaussian process of the same dimension as object. Any user function can be included with the only limitation that the two first parameters are:

• n: number of projections

• t: discretization points for functional problems

• m: number of columns for multivariate problems.

That functions must return a fdata or matrix object respectively.

The function allows user-defined contrasts. The list of contrast to be used for some of the factors in the formula. Each contrast matrix in the list has r rows, where r is the number of factor levels. The user can also request special predetermined contrasts, for example using the contr.helmert, contr.sum or contr.treatment functions.

The function returns (by default) the significance of the variables using the Bonferroni test and the False Discovery Rate test. Bootstrap procedure provides more precision

Value

An object with the following components:

• proj: The projection value of each point on the curves. Matrix with dimensions (RP x m), where RP is the number of projections and m is the number of points observed in each projection curve.

• mins: Minimum number for each random projection.

• result: p-value for each random projection.

• test.Bonf: Significance (TRUE or FALSE) for vector of random projections RP in columns and factors (and special contrasts) in rows.

• p.Bonf: p-value for vector of random projections RP in columns and factors (and special contrasts) in rows.

• test.fdr: False Discovery Rate (TRUE or FALSE) for vector of random projections RP in columns and factors (and special contrasts) in rows.

• p.fdr: p-value of False Discovery Rate for vector of random projections RP in columns and factors (and special contrasts) in rows.

• test.Boot: False Discovery Rate (TRUE or FALSE) for vector of random projections RP in columns and factors (and special contrasts) in rows.

• p.Boot: p-value of Bootstrap sample for vector of random projections RP in columns and factors (and special contrasts) in rows.

Note

anova.RPm deprecated.

If hetero=TRUE then all factors must be categorical.

Author(s)

Juan A. Cuesta-Albertos, Manuel Febrero-Bande, Manuel Oviedo de la Fuente
manuel.oviedo@udc.es

References

Cuesta-Albertos, J.A., Febrero-Bande, M. A simple multiway ANOVA for functional data. TEST 2010, DOI 10.1007/s11749-010-0185-3.

See Also

See Also as: fanova.onefactor

Examples

## Not run: 
# ex fanova.hetero
data(phoneme)
names(phoneme)
# A MV matrix obtained from functional data
data=as.data.frame(phoneme$learn$data[,c(1,seq(0,150,10)[-1])]) 
group=phoneme$classlearn
n=nrow(data)
group.rand=as.factor(sample(rep(1:3,len=n),n))
RP=c(2,5,15,30)

#ex 1: real factor and random factor
m03=data.frame(group,group.rand)
resul1=fanova.RPm(phoneme$learn,~group+group.rand,m03,RP=c(5,30))
summary(resul1)

#ex 2: real factor with special contrast
m0=data.frame(group)
cr5=contr.sum(5)   #each level vs last level
resul03c1=fanova.RPm(data,~group,m0,contrast=list(group=cr5))
summary(resul03c1)

#ex 3: random factor with special contrast. Same projs as ex 2.
m0=data.frame(group.rand)
zz=resul03c1$proj
cr3=contr.sum(3)   #each level vs last level
resul03c1=fanova.RPm(data,~group.rand,m0,contrast=list(group.rand=cr3),zproj=zz)
summary(resul03c1)

## End(Not run)
  

fda.usc internal functions

Description

Internal undocumentation functions for fda.usc package.

Usage

trace.matrix(x, na.rm = TRUE)

argvals.equi(tt)

## S3 method for class 'fdata'
fdata1 + fdata2

## S3 method for class 'fdata'
fdata1 - fdata2

## S3 method for class 'fdata'
fdata1 * fdata2

## S3 method for class 'fdata'
fdata1 / fdata2

## S3 method for class 'fdata'
fdataobj[i = TRUE, j = TRUE, drop = FALSE]

## S3 method for class 'fdata'
fdata1 != fdata2

## S3 method for class 'fdata'
fdata1 == fdata2

## S3 method for class 'fdata'
fdataobj ^ pot

## S3 method for class 'fdata'
dim(x)

ncol.fdata(x)

nrow.fdata(x)

## S3 method for class 'fdata'
length(x)

NROW.fdata(x)

NCOL.fdata(x)

rownames.fdata(x)

colnames.fdata(x)

## S3 method for class 'fdata'
c(...)

argvals(fdataobj)

rangeval(fdataobj)

## S3 method for class 'fdist'
fdataobj[i = TRUE, j = TRUE, drop = FALSE]

## S3 method for class 'fdata'
is.na(x)

## S3 method for class 'fdata'
anyNA(x, recursive = FALSE)

count.na.fdata(x)

unlist_fdata(x, recursive = TRUE, use.names = TRUE)

Arguments

Details

argvals.equi function returns TRUE if the argvals are equispaced and FALSE in othercase.

Note

In "Ops" functions "+.fdata", "-.fdata", "*.fdata" and "/.fdata": The lengths of the objects fdata1 and fdata2 may be different because operates recycled into minimum size as necessary.

References

Febrero-Bande, M., Oviedo de la Fuente, M. (2012). Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28. doi:10.18637/jss.v051.i04

Converts raw data or other functional data classes into fdata class.

Description

Create a functional data object of class fdata from (matrix, data.frame, numeric, integer, fd, fds, fts or sfts) class data.

Usage

fdata(mdata, argvals = NULL, rangeval = NULL, names = NULL, fdata2d = FALSE)

Arguments

Value

Return fdata class object with:

• "data": matrix of set cases with dimension (n x m), where n is the number of curves and m are the points observed in each curve

• "rangeval": the discretizations points values, if not provided: 1:m

• "rangeval": range of the discretizations points values, by default: range(argvals)

• "names": (optional) list with main an overall title, xlab title for x axis and ylab title for y axis.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Febrero-Bande, M., Oviedo de la Fuente, M. (2012). Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28. https://www.jstatsoft.org/v51/i04/

See Also

See Also as plot.fdata

Examples

## Not run: 
data(phoneme)
mlearn<-phoneme$learn[1:4,1:150]
# Center curves
fdata.c=fdata.cen(mlearn)$Xcen
par(mfrow=c(2,1))
plot(mlearn,type="l")
plot(fdata.c,type="l")

# Convert  from class fda to fdata
bsp1 <- create.bspline.basis(c(1,150),21)
fd1 <- Data2fd(1:150,y=t(mlearn$data),basisobj=bsp1)
fdataobj=fdata(fd1)

# Convert  from class fds, fts or sfts to fdata
#require(fds)
#a=fds(x = 1:20, y = Simulationdata$y, xname = "x", 
# yname = "Simulated value")
#b=fts(x = 15:49, y = Australiasmoothfertility$y, xname = "Age",
#    yname = "Fertility rate")
#c=sfts(ts(as.numeric(ElNino_ERSST_region_1and2$y), frequency = 12), xname = "Month",
#yname = "Sea surface temperature")
#class(a);class(b);class(c)
#fdataobj=fdata(b)

## End(Not run)

Bootstrap samples of a functional statistic

Description

provides bootstrap samples for functional data.

Usage

fdata.bootstrap(
  fdataobj,
  statistic = func.mean,
  alpha = 0.05,
  nb = 200,
  smo = 0,
  draw = FALSE,
  draw.control = NULL,
  ...
)

Arguments

Details

The fdata.bootstrap computes a confidence ball using bootstrap in the following way:

• Let X_1(t),\ldots,X_n(t) be the original data and T=T(X_1(t),\ldots,X_n(t)) be the sample statistic.

• Calculate the nb bootstrap resamples \left\{X_{1}^{*}(t),\cdots,X_n^*(t)\right\}, using the following scheme: X_i^*(t)=X_i(t)+Z(t), where Z(t) is normally distributed with mean 0 and covariance matrix \gamma\Sigma_x, where \Sigma_x is the covariance matrix of \left\{X_1(t),\ldots,X_n(t)\right\} and \gamma is the smoothing parameter.

• Let T^{*j}=T(X^{*j}_1(t),...,X^{*j}_n(t)) be the estimate using the j resample.

• Compute d(T,T^{*j}), j=1,\ldots,nb. Define the bootstrap confidence ball of level 1-\alpha as CB(\alpha)=X\in E such that d(T,X)\leq d_{\alpha} being d_{\alpha} the quantile (1-\alpha) of the distances between the bootstrap resamples and the sample estimate.

The fdata.bootstrap function allows us to define a statistic calculated on the nb resamples, control the degree of smoothing by smo argument and represent the confidence ball with level 1-\alpha as those resamples that fulfill the condition of belonging to CB(\alpha). The statistic used by default is the mean (func.mean) but also other depth-based functions can be used (see help(Descriptive)).

Value

• statistic: fdata class object with the statistic estimate from nb bootstrap samples.

• dband: Bootstrap estimate of (1-alpha)% distance.

• rep.dist: Distance from every replicate.

• resamples: fdata class object with the bootstrap resamples.

• fdataobj: fdata class object.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Cuevas A., Febrero-Bande, M. and Fraiman, R. (2007). Robust estimation and classification for functional data via projection-based depth notions. Computational Statistics 22, 3: 481-496.

Cuevas A., Febrero-Bande, M., Fraiman R. 2006. On the use of bootstrap for estimating functions with functional data. Computational Statistics and Data Analysis 51: 1063-1074.

Febrero-Bande, M., Oviedo de la Fuente, M. (2012). Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28. https://www.jstatsoft.org/v51/i04/

See Also

See Also as Descriptive

Examples

## Not run: 
data(tecator)
absorp<-tecator$absorp.fdata
# Time consuming
#Bootstrap for Trimmed Mean with depth mode
out.boot=fdata.bootstrap(absorp,statistic=func.trim.FM,nb=200,draw=TRUE)
names(out.boot)
#Bootstrap for Median with with depth mode
control=list("col"=c("grey","blue","cyan"),"lty"=c(2,1,1),"lwd"=c(1,3,1))
out.boot=fdata.bootstrap(absorp,statistic=func.med.mode,
draw=TRUE,draw.control=control)

## End(Not run)

Functional data centred (subtract the mean of each discretization point)

Description

The function fdata.cen centres the curves by subtracting the functional mean.

Usage

fdata.cen(fdataobj, meanX = func.mean(fdataobj))

Arguments

Value

Return:
two fdata class objects with:

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

See Also

See Also as fdata

Examples

## Not run: 
data(phoneme)
mlearn<-phoneme[["learn"]][13:15,]
fdata.c=fdata.cen(mlearn)$Xcen
par(mfrow=c(1,2))
plot(mlearn,type="l")
plot(fdata.c,type="l")

## End(Not run)

Computes the derivative of functional data object.

Description

Computes the derivative of functional data.

• If method ="bspline", "exponential", "fourier", "monomial" or "polynomial". fdata.deriv function creates a basis to represent the functional data. The functional data are converted to class fd using the Data2fd function and the basis indicated in the method. Finally, the function calculates the derivative of order nderiv of curves using deriv.fd function.
  

• If method="fmm", "periodic", "natural" or "monoH.FC" is used splinefun function.

• If method="diff", raw derivation is applied. Not recommended to use this method when the values are not equally spaced.
  

Usage

fdata.deriv(
  fdataobj,
  nderiv = 1,
  method = "bspline",
  class.out = "fdata",
  nbasis = NULL,
  ...
)

Arguments

Value

Returns the derivative of functional data of fd class if class.out="fd" or fdata class if class.out="fdata".

See Also

See also deriv.fd , splinefun and fdata

Examples

data(tecator)
absorp=tecator$absorp.fdata
tecator.fd1=fdata2fd(absorp)
tecator.fd2=fdata2fd(absorp,"fourier",9)
tecator.fd3=fdata2fd(absorp,"fourier",nbasis=9,nderiv=1)
#tecator.fd1;tecator.fd2;tecator.fd3
tecator.fdata1=fdata(tecator.fd1)
tecator.fdata2=fdata(tecator.fd2)
tecator.fdata3=fdata(tecator.fd3)
tecator.fdata4=fdata.deriv(absorp,nderiv=1,method="bspline",
class.out='fdata',nbasis=9)
tecator.fd4=fdata.deriv(tecator.fd3,nderiv=0,class.out='fd',nbasis=9)
plot(tecator.fdata4)
plot(fdata.deriv(absorp,nderiv=1,method="bspline",class.out='fd',nbasis=11))

fdata S3 Group Generic Functions

Description

fdata Group generic methods defined for four specified groups of functions, Math, Ops, Summary and Complex.

order.fdata and split.fdata: A wrapper for the order and split function for fdata object.

Usage

## S3 method for class 'fdata'
Math(x, ...)

## S3 method for class 'fdata'
Ops(e1, e2 = NULL)

## S3 method for class 'fdata'
Summary(..., na.rm = FALSE)

## S3 method for class 'fdata'
split(x, f, drop = FALSE, ...)

order.fdata(y, fdataobj, na.last = TRUE, decreasing = FALSE)

is.fdata(fdataobj)

Arguments

Details

In order.fdata the funcional data is ordered w.r.t the sample order of the values of vector.

split.fdata divides the data in the fdata object x into the groups defined by f.

Value

• split.fdata: The value returned from split is a list of fdata objects containing the values for the groups. The components of the list are named by the levels of f (after converting to a factor, or if already a factor and drop = TRUE, dropping unused levels).\

• order.fdata: returns the functional data fdataobj w.r.t. a permutation which rearranges its first argument into ascending or descending order.

Author(s)

Manuel Febrero Bande and Manuel Oviedo de la Fuente manuel.oviedo@udc.es

See Also

See Summary and Complex.

Examples

## Not run: 
data(tecator)
absor<-tecator$absorp.fdata
absor2<-fdata.deriv(absor,1)
absor<-absor2[1:5,1:4]
absor2<-absor2[1:5,1:4]
sum(absor)
round(absor,4)
log1<-log(absor)

fdataobj<-fdata(MontrealTemp)
fac<-factor(c(rep(1,len=17),rep(2,len=17)))
a1<-split(fdataobj,fac)
dim(a1[[1]]);dim(a1[[2]])

## End(Not run)

Compute fucntional coefficients from functional data represented in a base of functions

Description

Compute fucntional coefficients from functional data (fdata class object) represented in a basis (fixed of data-driven basis).

Usage

fdata2basis(fdataobj, basis, method = c("grid", "inprod"))

## S3 method for class 'basis.fdata'
summary(object, draw = TRUE, index = NULL, ...)

Arguments

Value

The fdata2basis function returns:

• coef: A matrix or two-dimensional array of coefficients.

• basis: Basis of fdata class evaluated on the same grid as fdataobj.

And summary function return:

• R: a matrix with a measure similar to R-sq for each curve aproximation (by row) and number of basis elements (by column).

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@usc.es

See Also

Inverse function: gridfdata. Alternative method: fdata2pc, fdata2pls.

Examples

## Not run: 
T <- 71
S <- 51
tj <- round(seq(0,1,len=T),3)
si <- round(seq(0,1,len=S),3)
beta1 <- outer(si,tj,function(si,tj){exp(-5*abs((tj-si)/5))})
nbasis.s =7
nbasis.t=11
base.s <- create.fourier.basis(c(0,1),nbasis=nbasis.s)
base.t <- create.fourier.basis(c(0,1),nbasis=nbasis.t)
y1 <- fdata(rbind(log(1+tj),1-5*(tj-0.5)^2),argvals=tj,rangeval=c(0,1))
aa <- fdata2basis(y1,base.t,method="inprod")
summary(aa)
plot(gridfdata(aa$coefs,aa$basis))
lines(y1,lwd=2,col=c(3,4),lty=2)

## End(Not run)

Converts fdata class object into fd class object

Description

Converts fdata class object into fd class object using Data2fd function.

Usage

fdata2fd(
  fdataobj,
  type.basis = NULL,
  nbasis = NULL,
  nderiv = 0,
  lambda = NULL,
  ...
)

Arguments

Value

Return an object of the fd class.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Febrero-Bande, M., Oviedo de la Fuente, M. (2012). Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28. https://www.jstatsoft.org/v51/i04/

See Also

See Also as fdata and Data2fd

Examples

## Not run: 
data(phoneme)
mlearn<-phoneme$learn[1]
fdata2=fdata2fd(mlearn)
class(mlearn)
class(fdata2)
fdata3=fdata2fd(mlearn,type.basis="fourier",nbasis=7)
plot(mlearn)
lines(fdata2,col=2)
lines(fdata3,col=3)
fdata5=fdata2fd(mlearn,nderiv=1)

## End(Not run)

Principal components for functional data

Description

Compute (penalized) principal components for functional data.

Usage

fdata2pc(fdataobj, ncomp = 2, norm = TRUE, lambda = 0, P = c(0, 0, 1), ...)

Arguments

Details

Smoothing is achieved by penalizing the integral of the square of the derivative of order m over rangeval:

• m = 0 penalizes the squared difference from 0 of the function

• m = 1 penalize the square of the slope or velocity

• m = 2 penalize the squared acceleration

• m = 3 penalize the squared rate of change of acceleration

Value

• d: The standard deviations of the functional principal components.

• rotation: Also known as loadings. A fdata class object whose rows contain the eigenvectors.

• x: Also known as scores. The value of the rotated functional data is returned.

• fdataobj.cen: The centered fdataobj object.

• mean: The functional mean of the fdataobj object.

• l: Vector of indices of principal components.

• C: The matched call.

• lambda: Amount of penalization.

• P: Penalty matrix.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@usc.es

References

Venables, W. N. and B. D. Ripley (2002). Modern Applied Statistics with S. Springer-Verlag.

N. Kraemer, A.-L. Boulsteix, and G. Tutz (2008). Penalized Partial Least Squares with Applications to B-Spline Transformations and Functional Data. Chemometrics and Intelligent Laboratory Systems, 94, 60 - 69. doi:10.1016/j.chemolab.2008.06.009

Febrero-Bande, M., Oviedo de la Fuente, M. (2012). Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28. https://www.jstatsoft.org/v51/i04/

See Also

See Also as svd and varimax.

Examples

 ## Not run: 
 n= 100;tt= seq(0,1,len=51)
 x0<-rproc2fdata(n,tt,sigma="wiener")
 x1<-rproc2fdata(n,tt,sigma=0.1)
 x<-x0*3+x1
 pc=fdata2pc(x,lambda=1)
 summary(pc)
 
## End(Not run)

Partial least squares components for functional data.

Description

Compute penalized partial least squares (PLS) components for functional data.

Usage

fdata2pls(fdataobj, y, ncomp = 2, lambda = 0, P = c(0, 0, 1), norm = TRUE, ...)

Arguments

Details

If norm=TRUE, computes the PLS by NIPALS algorithm and the Degrees of Freedom using the Krylov representation of PLS, see Kraemer and Sugiyama (2011).
If norm=FALSE, computes the PLS by Orthogonal Scores Algorithm and the Degrees of Freedom are the number of components ncomp, see Martens and Naes (1989).

Value

fdata2pls function return:

The fdata2pls function returns:

• df: Degree of freedom.

• rotation: fdata class object.

• x: The value of the rotated data (the centered data multiplied by the rotation matrix) is returned.

• fdataobj.cen: The centered fdataobj object.

• mean: Mean of fdataobj.

• l: Vector of indices of principal components.

• C: The matched call.

• lambda: Amount of penalization.

• P: Penalty matrix.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@usc.es

References

Kraemer, N., Sugiyama M. (2011). The Degrees of Freedom of Partial Least Squares Regression. Journal of the American Statistical Association. Volume 106, 697-705.

Febrero-Bande, M., Oviedo de la Fuente, M. (2012). Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28. https://www.jstatsoft.org/v51/i04/

Martens, H., Naes, T. (1989) Multivariate calibration. Chichester: Wiley.

See Also

Used in: fregre.pls, fregre.pls.cv. Alternative method: fdata2pc.

Examples

## Not run: 
n= 500;tt= seq(0,1,len=101)
x0<-rproc2fdata(n,tt,sigma="wiener")
x1<-rproc2fdata(n,tt,sigma=0.1)
x<-x0*3+x1
beta = tt*sin(2*pi*tt)^2
fbeta = fdata(beta,tt)
y<-inprod.fdata(x,fbeta)+rnorm(n,sd=0.1)
pls1=fdata2pls(x,y)
pls1$call
summary(pls1)
pls1$l
norm.fdata(pls1$rotation)

## End(Not run)

False Discorvery Rate (FDR)

Description

Compute the False Discovery Rate for a vector of p-values and alpha value.

Usage

FDR(pvalues = NULL, alpha = 0.95, dep = 1)

pvalue.FDR(pvalues = NULL, dep = 1)

Arguments

Details

FDR method is used for multiple hypothesis testing to correct problems of multiple contrasts.
If dep = 1, the tests are positively correlated, for example when many tests are the same contrast.
If dep < 1 the tests are negatively correlated.

Value

Return:

• out.FDR=TRUE: If there are significative differences.

• pv.FDR: p-value for False Discovery Rate test.

Author(s)

Febrero-Bande, M. and Oviedo de la Fuente, M.

References

Benjamini, Y., Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Annals of Statistics. 29 (4): 1165-1188. DOI:10.1214/aos/1013699998.

See Also

Function used in fanova.RPm

Examples

 p=seq(1:50)/1000
 FDR(p)
 pvalue.FDR(p)
 FDR(p,alpha=0.9999)
 FDR(p,alpha=0.9)
 FDR(p,alpha=0.9,dep=-1)

Tests for checking the equality of distributions between two functional populations.

Description

Three tests for the equality of distributions of two populations are provided. The null hypothesis is that the two populations are the same

Usage

XYRP.test(X.fdata, Y.fdata, nproj = 10, npc = 5, test = c("KS", "AD"))

MMD.test(
  X.fdata,
  Y.fdata,
  metric = "metric.lp",
  B = 1000,
  alpha = 0.95,
  kern = "RBF",
  ops.metric = list(lp = 2),
  draw = FALSE
)

MMDA.test(
  X.fdata,
  Y.fdata,
  metric = "metric.lp",
  B = 1000,
  alpha = 0.95,
  kern = "RBF",
  ops.metric = list(lp = 2),
  draw = FALSE
)

fEqDistrib.test(
  X.fdata,
  Y.fdata,
  metric = "metric.lp",
  method = c("Exch", "WildB"),
  B = 5000,
  ops.metric = list(lp = 2),
  iboot = FALSE
)

Arguments

Details

XYRP.test computes the p-values using random projections. Requires kSamples library. MMD.test computes Maximum Mean Discrepancy p-values using permutations (see Sejdinovic et al, (2013)) and MMDA.test does the same using an asymptotic approximation. fEqDistrib.test checks the equality of distributions using an embedding in a RKHS and two bootstrap approximations for calibration.

Value

A list with the following components by function:

• XYRP.test, FDR.pv: p-value using FDR, proj.pv: Matrix of p-values obtained for projections.

• MMD.test, MMDA.test: stat: Statistic, p.value: p-value, thresh: Threshold at level alpha.

• fEqDistrib.test, result: data.frame with columns Stat and p.value, Boot: data.frame with bootstrap replicas if iboot=TRUE.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.febrero@usc.es

References

Sejdinovic, D., Sriperumbudur, B., Gretton, A., Fukumizu, K. Equivalence of distance-based and RKHS-based statistics in Hypothesis Testing The Annals of Statistics, 2013. DOI 10.1214/13-AOS1140.

See Also

fmean.test.fdata, cov.test.fdata.

Examples

## Not run: 
tt=seq(0,1,len=51)
bet=0
mu1=fdata(10*tt*(1-tt)^(1+bet),tt)
mu2=fdata(10*tt^(1+bet)*(1-tt),tt) 
fsig=1
X=rproc2fdata(100,tt,mu1,sigma="vexponential",par.list=list(scale=0.2,theta=0.35))
Y=rproc2fdata(100,tt,mu2,sigma="vexponential",par.list=list(scale=0.2*fsig,theta=0.35))
fmean.test.fdata(X,Y,npc=-.98,draw=TRUE)
cov.test.fdata(X,Y,npc=5,draw=TRUE)
bet=0.1
mu1=fdata(10*tt*(1-tt)^(1+bet),tt)
mu2=fdata(10*tt^(1+bet)*(1-tt),tt) 
fsig=1.5
X=rproc2fdata(100,tt,mu1,sigma="vexponential",par.list=list(scale=0.2,theta=0.35))
Y=rproc2fdata(100,tt,mu2,sigma="vexponential",par.list=list(scale=0.2*fsig,theta=0.35))
fmean.test.fdata(X,Y,npc=-.98,draw=TRUE)
cov.test.fdata(X,Y,npc=5,draw=TRUE)
XYRP.test(X,Y,nproj=15)
MMD.test(X,Y,B=1000)
fEqDistrib.test(X,Y,B=1000)

## End(Not run)

Tests for checking the equality of means and/or covariance between two populations under gaussianity.

Description

Two tests for the equality of means and covariances of two populations are provided. Both tests are constructed under gaussianity following Horvath & Kokoszka, 2012, Chapter 5.

Usage

fmean.test.fdata(
  X.fdata,
  Y.fdata,
  method = c("X2", "Boot"),
  npc = 5,
  alpha = 0.95,
  B = 1000,
  draw = FALSE
)

cov.test.fdata(
  X.fdata,
  Y.fdata,
  method = c("X2", "Boot"),
  npc = 5,
  alpha = 0.95,
  B = 1000,
  draw = FALSE
)

Arguments

Details

fmean.test.fdata computes the test for equality of means. cov.test.fdata computes the test for equality of covariance operators. Both tests have asymptotic distributions under the null related with chi-square distribution. Also, a parametric bootstrap procedure is implemented in both cases.

Value

Return a list with:

• stat: Value of the statistic.

• pvalue: P-values for the test.

• vcrit: Critical cutoff for rejecting the null hypothesis.

• p: Degrees of freedom for X2 statistic.

• B: Number of bootstrap replicas.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.febrero@usc.es

References

Inference for Functional Data with Applications. Horvath, L and Kokoszka, P. (2012). Springer.

See Also

See Also as fanova.RPm, fanova.onefactor.

Examples

## Not run: 
tt=seq(0,1,len=51)
bet=0
mu1=fdata(10*tt*(1-tt)^(1+bet),tt)
mu2=fdata(10*tt^(1+bet)*(1-tt),tt) 
fsig=1
X=rproc2fdata(100,tt,mu1,sigma="vexponential",par.list=list(scale=0.2,theta=0.35))
Y=rproc2fdata(100,tt,mu2,sigma="vexponential",par.list=list(scale=0.2*fsig,theta=0.35))
fmean.test.fdata(X,Y,npc=-.98,draw=TRUE)
cov.test.fdata(X,Y,npc=5,draw=TRUE)
bet=0.1
mu1=fdata(10*tt*(1-tt)^(1+bet),tt)
mu2=fdata(10*tt^(1+bet)*(1-tt),tt) 
fsig=1.5
X=rproc2fdata(100,tt,mu1,sigma="vexponential",par.list=list(scale=0.2,theta=0.35))
Y=rproc2fdata(100,tt,mu2,sigma="vexponential",par.list=list(scale=0.2*fsig,theta=0.35))
fmean.test.fdata(X,Y,npc=-.98,draw=TRUE)
cov.test.fdata(X,Y,npc=5,draw=TRUE)

## End(Not run)

F-test for the Functional Linear Model with scalar response

Description

The function flm.Ftest tests the null hypothesis of no interaction between a functional covariate and a scalar response inside the Functional Linear Model (FLM): Y=\big<X,\beta\big>+\epsilon. The null hypothesis is H_0:\,\beta=0 and the alternative is H_1:\,\beta\neq 0. The null hypothesis is tested by a functional extension of the classical F-test (see Details).

Usage

Ftest.statistic(X.fdata, Y)

flm.Ftest(X.fdata, Y, B = 5000, verbose = TRUE)

Arguments

Details

The Functional Linear Model with scalar response (FLM), is defined as Y=\big<X,\beta\big>+\epsilon, for a functional process X such that E[X(t)]=0, E[X(t)\epsilon]=0 for all t and for a scalar variable Y such that E[Y]=0. The functional F-test is defined as

T_n=\bigg\|\frac{1}{n}\sum_{i=1}^n (X_i-\bar X)(Y_i-\bar Y)\bigg\|,

where \bar X is the functional mean of X, \bar Y is the ordinary mean of Y and \|\cdot\| is the L^2 functional norm. The statistic is computed with the function Ftest.statistic. The distribution of the test statistic is approximated by a wild bootstrap resampling on the residuals, using the golden section bootstrap.

Value

The value for Ftest.statistic is simply the F-test statistic. The value for flm.Ftest is an object with class "htest" whose underlying structure is a list containing the following components:

• statistic: The value of the F-test statistic.

• boot.statistics: A vector of length B with the values of the bootstrap F-test statistics.

• p.value: The p-value of the test.

• method: The character string "Functional Linear Model F-test".

• B: The number of bootstrap replicates used.

• data.name: The character string "Y=<X,0>+e".

Note

No NA's are allowed neither in the functional covariate nor in the scalar response.

Author(s)

Eduardo Garcia-Portugues. Please, report bugs and suggestions to eduardo.garcia.portugues@uc3m.es

References

Garcia-Portugues, E., Gonzalez-Manteiga, W. and Febrero-Bande, M. (2014). A goodness–of–fit test for the functional linear model with scalar response. Journal of Computational and Graphical Statistics, 23(3), 761-778. doi:10.1080/10618600.2013.812519

Gonzalez-Manteiga, W., Gonzalez-Rodriguez, G., Martinez-Calvo, A. and Garcia-Portugues, E. Bootstrap independence test for functional linear models. arXiv:1210.1072. https://arxiv.org/abs/1210.1072

See Also

rwild, flm.test, dfv.test

Examples

## Not run: 
## Simulated example ##
X=rproc2fdata(n=50,t=seq(0,1,l=101),sigma="OU")
beta0=fdata(mdata=rep(0,length=101)+rnorm(101,sd=0.05),
argvals=seq(0,1,l=101),rangeval=c(0,1))
beta1=fdata(mdata=cos(2*pi*seq(0,1,l=101))-(seq(0,1,l=101)-0.5)^2+
rnorm(101,sd=0.05),argvals=seq(0,1,l=101),rangeval=c(0,1))

# Null hypothesis holds
Y0=drop(inprod.fdata(X,beta0)+rnorm(50,sd=0.1))
# Null hypothesis does not hold
Y1=drop(inprod.fdata(X,beta1)+rnorm(50,sd=0.1))

# Do not reject H0
flm.Ftest(X,Y0,B=100)
flm.Ftest(X,Y0,B=5000)

# Reject H0
flm.Ftest(X,Y1,B=100)
flm.Ftest(X,Y1,B=5000)

## End(Not run)

Goodness-of-fit test for the Functional Linear Model with scalar response

Description

The function flm.test tests the composite null hypothesis of a Functional Linear Model with scalar response (FLM),

H_0:\,Y=\big<X,\beta\big>+\epsilon,

versus a general alternative. If \beta=\beta_0 is provided, then the simple hypothesis H_0:\,Y=\big<X,\beta_0\big>+\epsilon is tested. The testing of the null hypothesis is done by a Projected Cramer-von Mises statistic (see Details).

Usage

flm.test(
  X.fdata,
  Y,
  beta0.fdata = NULL,
  B = 5000,
  est.method = "pls",
  p = NULL,
  type.basis = "bspline",
  verbose = TRUE,
  plot.it = TRUE,
  B.plot = 100,
  G = 200,
  ...
)

Arguments

Details

The Functional Linear Model with scalar response (FLM), is defined as Y=\big<X,\beta\big>+\epsilon, for a functional process X such that E[X(t)]=0, E[X(t)\epsilon]=0 for all t and for a scalar variable Y such that E[Y]=0. Then, the test assumes that Y and X.fdata are centred and will automatically center them. So, bear in mind that when you apply the test for Y and X.fdata, actually, you are applying it to Y-mean(Y) and fdata.cen(X.fdata)$Xcen. The test statistic corresponds to the Cramer-von Mises norm of the Residual Marked empirical Process based on Projections R_n(u,\gamma) defined in Garcia-Portugues et al. (2014). The expression of this process in a p-truncated basis of the space L^2[0,T] leads to the p-multivariate process R_{n,p}\big(u,\gamma^{(p)}\big), whose Cramer-von Mises norm is computed. The choice of an appropriate p to represent the functional process X, in case that is not provided, is done via the estimation of \beta for the composite hypothesis. For the simple hypothesis, as no estimation of \beta is done, the choice of p depends only on the functional process X. As the result of the test may change for different p's, we recommend to use an automatic criterion to select p instead of provide a fixed one. The distribution of the test statistic is approximated by a wild bootstrap resampling on the residuals, using the golden section bootstrap. Finally, the graph shown if plot.it=TRUE represents the observed trajectory, and the bootstrap trajectories under the null, of the process RMPP integrated on the projections:

R_n(u)\approx\frac{1}{G}\sum_{g=1}^G R_n(u,\gamma_g),

where \gamma_g are simulated as Gaussians processes. This gives a graphical idea of how distant is the observed trajectory from the null hypothesis.

Value

An object with class "htest" whose underlying structure is a list containing the following components:

• statistic: The value of the test statistic.

• boot.statistics: A vector of length B with the values of the bootstrap test statistics.

• p.value: The p-value of the test.

• method: The method used.

• B: The number of bootstrap replicates used.

• type.basis: The type of basis used.

• beta.est: The estimated functional parameter \beta in the composite hypothesis. For the simple hypothesis, the given beta0.fdata.

• p: The number of basis elements passed or automatically chosen.

• ord: The optimal order for PC and PLS given by fregre.pc.cv and fregre.pls.cv. For other methods, it is set to 1:p.

• data.name: The character string "Y=<X,b>+e".

Note

No NA's are allowed neither in the functional covariate nor in the scalar response.

Author(s)

Eduardo Garcia-Portugues. Please, report bugs and suggestions to edgarcia@est-econ.uc3m.es

References

Escanciano, J. C. (2006). A consistent diagnostic test for regression models using projections. Econometric Theory, 22, 1030-1051. doi:10.1017/S0266466606060506

Garcia-Portugues, E., Gonzalez-Manteiga, W. and Febrero-Bande, M. (2014). A goodness–of–fit test for the functional linear model with scalar response. Journal of Computational and Graphical Statistics, 23(3), 761-778. doi:10.1080/10618600.2013.812519

See Also

Adot, PCvM.statistic, rwild, flm.Ftest, dfv.test, fregre.pc, fregre.pls, fregre.basis, fregre.pc.cv, fregre.pls.cv, fregre.basis.cv, optim.basis, create.basis

Examples

# Simulated example #
X=rproc2fdata(n=100,t=seq(0,1,l=101),sigma="OU")
beta0=fdata(mdata=cos(2*pi*seq(0,1,l=101))-(seq(0,1,l=101)-0.5)^2+
            rnorm(101,sd=0.05),argvals=seq(0,1,l=101),rangeval=c(0,1))
Y=inprod.fdata(X,beta0)+rnorm(100,sd=0.1)

dev.new(width=21,height=7)
par(mfrow=c(1,3))
plot(X,main="X")
plot(beta0,main="beta0")
plot(density(Y),main="Density of Y",xlab="Y",ylab="Density")
rug(Y)

## Not run: 
# Composite hypothesis: do not reject FLM
pcvm.sim=flm.test(X,Y,B=50,B.plot=50,G=100,plot.it=TRUE)
pcvm.sim
flm.test(X,Y,B=5000)
 
# Estimated beta
dev.new()
plot(pcvm.sim$beta.est)

# Simple hypothesis: do not reject beta=beta0
flm.test(X,Y,beta0.fdata=beta0,B=50,B.plot=50,G=100)
flm.test(X,Y,beta0.fdata=beta0,B=5000) 

# AEMET dataset #
data(aemet)
# Remove the 5\
dev.new()
res.FM=depth.FM(aemet$temp,draw=TRUE)
qu=quantile(res.FM$dep,prob=0.05)
l=which(res.FM$dep<=qu)
lines(aemet$temp[l],col=3)
aemet$df$name[l]

# Data without outliers 
wind.speed=apply(aemet$wind.speed$data,1,mean)[-l]
temp=aemet$temp[-l]
# Exploratory analysis: accept the FLM
pcvm.aemet=flm.test(temp,wind.speed,est.method="pls",B=100,B.plot=50,G=100)
pcvm.aemet

# Estimated beta
dev.new()
plot(pcvm.aemet$beta.est,lwd=2,col=2)
# B=5000 for more precision on calibration of the test: also accept the FLM
flm.test(temp,wind.speed,est.method="pls",B=5000) 

# Simple hypothesis: rejection of beta0=0? Limiting p-value...
dat=rep(0,length(temp$argvals))
flm.test(temp,wind.speed, beta0.fdata=fdata(mdata=dat,argvals=temp$argvals,
                                            rangeval=temp$rangeval),B=100)
flm.test(temp,wind.speed, beta0.fdata=fdata(mdata=dat,argvals=temp$argvals,
                                            rangeval=temp$rangeval),B=5000) 
                                            
# Tecator dataset #
data(tecator)
names(tecator)
absorp=tecator$absorp.fdata
ind=1:129 # or ind=1:215
x=absorp[ind,]
y=tecator$y$Fat[ind]
tt=absorp[["argvals"]]

# Exploratory analysis for composite hypothesis with automatic choose of p
pcvm.tecat=flm.test(x,y,B=100,B.plot=50,G=100)
pcvm.tecat

# B=5000 for more precision on calibration of the test: also reject the FLM
flm.test(x,y,B=5000) 

# Distribution of the PCvM statistic
plot(density(pcvm.tecat$boot.statistics),lwd=2,xlim=c(0,10),
              main="PCvM distribution", xlab="PCvM*",ylab="Density")
rug(pcvm.tecat$boot.statistics)
abline(v=pcvm.tecat$statistic,col=2,lwd=2)
legend("top",legend=c("PCvM observed"),lwd=2,col=2)

# Simple hypothesis: fixed p
dat=rep(0,length(x$argvals))
flm.test(x,y,beta0.fdata=fdata(mdata=dat,argvals=x$argvals,
                               rangeval=x$rangeval),B=100,p=11)
                               
# Simple hypothesis, automatic choose of p
flm.test(x,y,beta0.fdata=fdata(mdata=dat,argvals=x$argvals,
                               rangeval=x$rangeval),B=100)
flm.test(x,y,beta0.fdata=fdata(mdata=dat,argvals=x$argvals,
                               rangeval=x$rangeval),B=5000)

## End(Not run)

Functional Regression with scalar response using basis representation.

Description

Computes functional regression between functional explanatory variable X(t) and scalar response Y using basis representation.

Usage

fregre.basis(
  fdataobj,
  y,
  basis.x = NULL,
  basis.b = NULL,
  lambda = 0,
  Lfdobj = vec2Lfd(c(0, 0), rtt),
  weights = rep(1, n),
  ...
)

Arguments

Details

Y=\big<X,\beta\big>+\epsilon=\int_{T}{X(t)\beta(t)dt+\epsilon}

where \big< \cdot , \cdot \big> denotes the inner product on L_2 and \epsilon are random errors with mean zero, finite variance \sigma^2 and E[X(t)\epsilon]=0.

The function uses the basis representation proposed by Ramsay and Silverman (2005) to model the relationship between the scalar response and the functional covariate by basis representation of the observed functional data X(t)\approx\sum_{k=1}^{k_{n1}} c_k \xi_k(t) and the unknown functional parameter \beta(t)\approx\sum_{k=1}^{k_{n2}} b_k \phi_k(t).

The functional linear models estimated by the expression:

\hat{y}= \big< X,\hat{\beta} \big> = C^{T}\psi(t)\phi^{T}(t)\hat{b}=\tilde{X}\hat{b}

where \tilde{X}(t)=C^{T}\psi(t)\phi^{T}(t), and \hat{b}=(\tilde{X}^{T}\tilde{X})^{-1}\tilde{X}^{T}y and so, \hat{y}=\tilde{X}\hat{b}=\tilde{X}(\tilde{X}^{T}\tilde{X})^{-1}\tilde{X}^{T}y=Hy where H is the hat matrix with degrees of freedom: df=tr(H).

If \lambda>0 then fregre.basis incorporates a roughness penalty:
\hat{y}=\tilde{X}\hat{b}=\tilde{X}(\tilde{X}^{T}\tilde{X}+\lambda R_0)^{-1}\tilde{X}^{T}y= H_{\lambda}y where R_0 is the penalty matrix.

This function allows covariates of class fdata, matrix, data.frame or directly covariates of class fd. The function also gives default values to arguments basis.x and basis.b for representation on the basis of functional data X(t) and the functional parameter \beta(t), respectively.

If basis=NULL creates the bspline basis by create.bspline.basis.
If the functional covariate fdataobj is a matrix or data.frame, it creates an object of class "fdata" with default attributes, see fdata.
If basis.x$type=``fourier'' and basis.b$type=``fourier'', the basis are orthonormal and the function decreases the number of fourier basis elements on the min(k_{n1},k_{n2}), where k_{n1} and k_{n2} are the number of basis element of basis.x and basis.b respectively.

Value

Return:

• call: The matched call.

• coefficients: A named vector of coefficients.

• residuals: y minus fitted values.

• fitted.values: Estimated scalar response.

• beta.est: Estimated beta parameter of class fd.

• weights: (only for weighted fits) the specified weights.

• df.residual: The residual degrees of freedom.

• r2: Coefficient of determination.

• sr2: Residual variance.

• Vp: Estimated covariance matrix for the parameters.

• H: Hat matrix.

• y: Response.

• fdataobj: Functional explanatory data of class fdata.

• a.est: Intercept parameter estimated.

• x.fd: Centered functional explanatory data of class fd.

• basis.b: Basis used for beta parameter estimation.

• lambda.opt: A roughness penalty.

• Lfdobj: Order of a derivative or a linear differential operator.

• P: Penalty matrix.

• lm: Return lm object.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Ramsay, James O., and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York.

Febrero-Bande, M., Oviedo de la Fuente, M. (2012). Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28. https://www.jstatsoft.org/v51/i04/

See Also

See Also as: fregre.basis.cv, summary.fregre.fd and predict.fregre.fd.
Alternative method: fregre.pc and fregre.np.

Examples

## Not run: 
# fregre.basis
data(tecator)
names(tecator)
absorp=tecator$absorp.fdata
ind=1:129
x=absorp[ind,]
y=tecator$y$Fat[ind]
tt=absorp[["argvals"]]
res1=fregre.basis(x,y)
summary(res1)
basis1=create.bspline.basis(rangeval=range(tt),nbasis=19)
basis2=create.bspline.basis(rangeval=range(tt),nbasis=9)
res5=fregre.basis(x,y,basis1,basis2)
summary(res5)
x.d2=fdata.deriv(x,nbasis=19,nderiv=1,method="bspline",class.out="fdata")
res7=fregre.basis(x.d2,y,basis1,basis2)
summary(res7)

## End(Not run)

Cross-validation Functional Regression with scalar response using basis representation.

Description

Computes functional regression between functional explanatory variables and scalar response using basis representation.

Usage

fregre.basis.cv(
  fdataobj,
  y,
  basis.x = NULL,
  basis.b = NULL,
  type.basis = NULL,
  lambda = 0,
  Lfdobj = vec2Lfd(c(0, 0), rtt),
  type.CV = GCV.S,
  par.CV = list(trim = 0),
  weights = rep(1, n),
  verbose = FALSE,
  ...
)

Arguments

Details

The function fregre.basis.cv() uses validation criterion defined by argument type.CV to estimate the number of basis elements and/or the penalized parameter (lambda) that best predicts the response.

If basis = NULL creates bspline basis.

If the functional covariate fdataobj is in a format raw data, such as matrix or data.frame, creates an object of class fdata with default attributes, see fdata.

If basis.x is a vector of number of basis elements and basis.b=NULL, the function force the same number of elements in the basis of x and beta.

If basis.x$type=``fourier'' and basis.b$type=``fourier'', the function decreases the number of fourier basis elements on the min(k_{n1},k_{n2}), where k_{n1} and k_{n2} are the number of basis element of basis.x and basis.b respectively.

Value

Return:

• fregre.basis: Fitted regression object by the best parameters (basis elements for data and beta and lambda penalty).

• basis.x.opt: Basis used for functional explanatory data estimation fdata.

• basis.b.opt: Basis used for functional beta parameter estimation.

• lambda.opt: lambda value that minimizes CV or GCV method.

• gcv.opt: Minimum value of CV or GCV method.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Ramsay, James O. and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York.

Febrero-Bande, M., Oviedo de la Fuente, M. (2012). Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28. https://www.jstatsoft.org/v51/i04/

See Also

See Also as: fregre.basis, summary.fregre.fd and predict.fregre.fd .
Alternative method: fregre.pc.cv and fregre.np.cv.

Examples

## Not run: 
data(tecator)
x<-tecator$absorp.fdata[1:129]
y=tecator$y$Fat[1:129]
b1<-c(15,21,31)
b2<-c(7,9)
res1=fregre.basis.cv(x,y,basis.x=b1)
res2=fregre.basis.cv(x,y,basis.x=b1,basis.b=b2)
res1$gcv
res2$gcv
l=2^(-4:10)
res3=fregre.basis.cv(x,y,basis.b=b1,type.basis="fourier",
lambda=l,type.CV=GCV.S,par.CV=list(trim=0.15))
res3$gcv

## End(Not run)

Functional Regression with functional response using basis representation.

Description

Computes functional regression between functional explanatory variable X(s) and functional response Y(t) using basis representation.

Usage

fregre.basis.fr(
  x,
  y,
  basis.s = NULL,
  basis.t = NULL,
  lambda.s = 0,
  lambda.t = 0,
  Lfdobj.s = vec2Lfd(c(0, 0), range.s),
  Lfdobj.t = vec2Lfd(c(0, 0), range.t),
  weights = NULL,
  ...
)

Arguments

Details

Y(t)=\alpha(t)+\int_{T}{X(s)\beta(s,t)ds+\epsilon(t)}

where \alpha(t) is the intercept function, \beta(s,t) is the bivariate resgression function and \epsilon(t) are the error term with mean zero.

The function is a wrapped of linmod function proposed by Ramsay and Silverman (2005) to model the relationship between the functional response Y(t) and the functional covariate X(t) by basis representation of both.

The unknown bivariate functional parameter \beta(s,t) can be expressed as a double expansion in terms of K basis function \nu_k and L basis functions \theta_l,

\beta(s,t)=\sum_{k=1}^{K}\sum_{l=1}^{L} b_{kl} \nu_{k}(s)\theta_{l}(t)=\nu(s)^{\top}\bold{B}\theta(t)

Then, the model can be re–written in a matrix version as,

Y(t)=\alpha(t)+\int_{T}{X(s)\nu(s)^{\top}\bold{B}\theta(t)ds+\epsilon(t)}=\alpha(t)+\bold{XB}\theta(t)+\epsilon(t)

where \bold{X}=\int X(s)\nu^{\top}(t)ds

This function allows objects of class fdata or directly covariates of class fd. If x is a fdata class, basis.s is also the basis used to represent x as fd class object. If y is a fdata class, basis.t is also the basis used to represent y as fd class object. The function also gives default values to arguments basis.s and basis.t for construct the bifd class object used in the estimation of \beta(s,t). If basis.s=NULL or basis.t=NULL the function creates a bspline basis by create.bspline.basis.

fregre.basis.fr incorporates a roughness penalty using an appropiate linear differential operator; lambda.s, Lfdobj.s for penalization of \beta's variations with respect to s and
lambda.t, Lfdobj.t for penalization of \beta's variations with respect to t.

Value

Return:

• call: The matched call.

• a.est: Intercept parameter estimated.

• coefficients: The matrix of the coefficients.

• beta.est: A bivariate functional data object of class bifd with the estimated parameters of \beta(s,t).

• fitted.values: Estimated response.

• residuals: y minus fitted values.

• y: Functional response.

• x: Functional explanatory data.

• lambda.s: A roughness penalty with respect to s.

• lambda.t: A roughness penalty with respect to t.

• Lfdobj.s: A linear differential operator with respect to s.

• Lfdobj.t: A linear differential operator with respect to t.

• weights: Weights.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Ramsay, James O., and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York.

See Also

See Also as: predict.fregre.fr. Alternative method: linmod.

Examples

## Not run: 
rtt<-c(0, 365)
basis.alpha  <- create.constant.basis(rtt)
basisx  <- create.bspline.basis(rtt,11)
basisy  <- create.bspline.basis(rtt,11)
basiss  <- create.bspline.basis(rtt,7)
basist  <- create.bspline.basis(rtt,9)

# fd class
dayfd<-Data2fd(day.5,CanadianWeather$dailyAv,basisx)
tempfd<-dayfd[,1]
log10precfd<-dayfd[,3]
res1 <-  fregre.basis.fr(tempfd, log10precfd,
basis.s=basiss,basis.t=basist)

# fdata class
tt<-1:365
tempfdata<-fdata(t(CanadianWeather$dailyAv[,,1]),tt,rtt)
log10precfdata<-fdata(t(CanadianWeather$dailyAv[,,3]),tt,rtt)
res2<-fregre.basis.fr(tempfdata,log10precfdata,
basis.s=basiss,basis.t=basist)

# penalization
Lfdobjt <- Lfdobjs <- vec2Lfd(c(0,0), rtt)
Lfdobjt <- vec2Lfd(c(0,0), rtt)
lambdat<-lambdas <- 100
res1.pen <- fregre.basis.fr(tempfdata,log10precfdata,basis.s=basiss,
basis.t=basist,lambda.s=lambdas,lambda.t=lambdat,
Lfdobj.s=Lfdobjs,Lfdobj.t=Lfdobjt)

res2.pen <- fregre.basis.fr(tempfd, log10precfd,
basis.s=basiss,basis.t=basist,lambda.s=lambdas,
lambda.t=lambdat,Lfdobj.s=Lfdobjs,Lfdobj.t=Lfdobjt)

plot(log10precfd,col=1)
lines(res1$fitted.values,col=2)
plot(res1$residuals)
plot(res1$beta.est,tt,tt)
plot(res1$beta.est,tt,tt,type="persp",theta=45,phi=30)

## End(Not run)

Bootstrap regression

Description

Estimate the beta parameter by wild or smoothed bootstrap procedure

Usage

fregre.bootstrap(
  model,
  nb = 500,
  wild = TRUE,
  type.wild = "golden",
  newX = NULL,
  smo = 0.1,
  smoX = 0.05,
  alpha = 0.95,
  kmax.fix = FALSE,
  draw = TRUE,
  ...
)

Arguments

Details

Estimate the beta parameter by wild or smoothed bootstrap procedure using principal components representation fregre.pc, Partial least squares components (PLS) representation fregre.pls or basis representation fregre.basis.
If a new curves are in newX argument the bootstrap method estimates the response using the bootstrap resamples.

If the model exhibits heteroskedasticity, the use of wild bootstrap procedure is recommended (by default).

Value

Return:

• model: fregre.pc, fregre.pls or fregre.basis object.

• beta.boot: Functional beta estimated by the nb bootstrap regressions.

• norm.boot: Norm of differences between the nboot betas estimated by bootstrap and beta estimated by the regression model.

• coefs.boot: Matrix with the bootstrap estimated basis coefficients.

• kn.boot: Vector or list of length nb with index of the basis, PC or PLS factors selected in each bootstrap regression.

• y.pred: Predicted response values using newX covariates.

• y.boot: Matrix of bootstrap predicted response values using newX covariates.

• newX: A fdata class containing the values of the model covariates at which predictions are required (only for smoothed bootstrap).

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Febrero-Bande, M., Galeano, P. and Gonzalez-Manteiga, W. (2010). Measures of influence for the functional linear model with scalar response. Journal of Multivariate Analysis 101, 327-339.

Febrero-Bande, M., Oviedo de la Fuente, M. (2012). Statistical Computing in Functional Data Analysis: The R Package fda.usc. Journal of Statistical Software, 51(4), 1-28. https://www.jstatsoft.org/v51/i04/

See Also

See Also as: fregre.pc, fregre.pls, fregre.basis, .

Examples

## Not run:  
data(tecator)
iest<-1:165
x=tecator$absorp.fdata[iest]
y=tecator$y$Fat[iest]
nb<-25  ## Time-consuming
res.pc=fregre.pc(x,y,1:6)
# Fix the compontents used in the each regression
res.boot1=fregre.bootstrap(res.pc,nb=nb,wild=FALSE,kmax.fix=TRUE)
# Select the "best" compontents used in the each regression
res.boot2=fregre.bootstrap(res.pc,nb=nb,wild=FALSE,kmax.fix=FALSE) 
res.boot3=fregre.bootstrap(res.pc,nb=nb,wild=FALSE,kmax.fix=10) 
## predicted responses and bootstrap confidence interval
newx=tecator$absorp.fdata[-iest]
res.boot4=fregre.bootstrap(res.pc,nb=nb,wild=FALSE,newX=newx,draw=TRUE)

## End(Not run)

Fitting Functional Generalized Kernel Additive Models.

Description

Computes functional regression between functional explanatory variables (X^{1}(t_1),...,X^{q}(t_q)) and scalar response Y using backfitting algorithm.

Usage

fregre.gkam(
  formula,
  family = gaussian(),
  data,
  weights = rep(1, nobs),
  par.metric = NULL,
  par.np = NULL,
  offset = NULL,
  control = list(maxit = 100, epsilon = 0.001, trace = FALSE, inverse = "solve"),
  ...
)

Arguments

Details

The smooth functions f(.) are estimated nonparametrically using a iterative local scoring algorithm by applying Nadaraya-Watson weighted kernel smoothers using fregre.np.cv in each step, see Febrero-Bande and Gonzalez-Manteiga (2011) for more details.
Consider the fitted response \hat{Y}=g^{-1}(H_{Q}y), where H_{Q} is the weighted hat matrix.
Opsomer and Ruppert (1997) solves a system of equations for fit the unknowns f(\cdot) computing the additive smoother matrix H_k such that \hat{f}_k (X^k)=H_{k}Y and H_Q=H_1+,\cdots,+H_q. The additive model is fitted as follows:

\hat{Y}=g^{-1}\Big(\sum_i^q \hat{f_i}(X_i)\Big)

Value

• result: List of non-parametric estimation by covariate.

• fitted.values: Estimated scalar response.

• residuals: y minus fitted values.

• effects: The residual degrees of freedom.

• alpha: Hat matrix.

• family: Coefficient of determination.

• linear.predictors: Residual variance.

• deviance: Scalar response.

• aic: Functional explanatory data.

• null.deviance: Non functional explanatory data.

• iter: Distance matrix between curves.

• w: Beta coefficient estimated.

• eqrank: List that containing the variables in the model.

• prior.weights: Asymmetric kernel used.

• y: Scalar response.

• H: Hat matrix, see Opsomer and Ruppert (1997) for more details.

• converged: Conv.

Author(s)

Febrero-Bande, M. and Oviedo de la Fuente, M.

References

Febrero-Bande M. and Gonzalez-Manteiga W. (2012). Generalized Additive Models for Functional Data. TEST. Springer-Velag. doi:10.1007/s11749-012-0308-0

Opsomer J.D. and Ruppert D.(1997). Fitting a bivariate additive model by local polynomial regression.Annals of Statistics, 25, 186-211.

See Also

See Also as: fregre.gsam, fregre.glm and fregre.np.cv

Examples

## Not run: 
data(tecator)
ab=tecator$absorp.fdata[1:100]
ab2=fdata.deriv(ab,2)
yfat=tecator$y[1:100,"Fat"]

# Example 1: # Changing the argument par.np and family
yfat.cat=ifelse(yfat<15,0,1)
xlist=list("df"=data.frame(yfat.cat),"ab"=ab,"ab2"=ab2)
f2<-yfat.cat~ab+ab2

par.NP<-list("ab"=list(Ker=AKer.norm,type.S="S.NW"),
"ab2"=list(Ker=AKer.norm,type.S="S.NW"))
res2=fregre.gkam(f2,family=binomial(),data=xlist,
par.np=par.NP)
res2

# Example 2: Changing the argument par.metric and family link
par.metric=list("ab"=list(metric=semimetric.deriv,nderiv=2,nbasis=15),
"ab2"=list("metric"=semimetric.basis))
res3=fregre.gkam(f2,family=binomial("probit"),data=xlist,
par.metric=par.metric,control=list(maxit=2,trace=FALSE))
summary(res3)

# Example 3: Gaussian family (by default)
# Only 1 iteration (by default maxit=100)
xlist=list("df"=data.frame(yfat),"ab"=ab,"ab2"=ab2)
f<-yfat~ab+ab2
res=fregre.gkam(f,data=xlist,control=list(maxit=1,trace=FALSE))
res

## End(Not run)

Fitting Functional Generalized Linear Models

Description

Computes functional generalized linear model between functional covariate X^j(t) (and non functional covariate Z^j) and scalar response Y using basis representation.

Usage

fregre.glm(
  formula,
  family = gaussian(),
  data,
  basis.x = NULL,
  basis.b = NULL,
  subset = NULL,
  weights = NULL,
  ...
)

Arguments

Details

This function is an extension of the linear regression models: fregre.lm where the E[Y|X,Z] is related to the linear prediction \eta via a link function g(.).

E[Y|X,Z]=\eta=g^{-1}(\alpha+\sum_{j=1}^{p}\beta_{j}Z^{j}+\sum_{k=1}^{q}\frac{1}{\sqrt{T_k}}\int_{T_k}{X^{k}(t)\beta_{k}(t)dt})

where Z=\left[ Z^1,\cdots,Z^p \right] are the non functional covariates and X(t)=\left[ X^{1}(t_1),\cdots,X^{q}(t_q) \right] are the functional ones.

The first item in the data list is called "df" and is a data frame with the response and non functional explanatory variables, as glm.

Functional covariates of class fdata or fd are introduced in the following items in the data list.
basis.x is a list of basis for represent each functional covariate. The basis object can be created by the function: create.pc.basis, pca.fd create.pc.basis, create.fdata.basis o create.basis.
basis.b is a list of basis for represent each \beta(t) parameter. If basis.x is a list of functional principal components basis (see create.pc.basis or pca.fd) the argument basis.b is ignored.

represent beta lower than the number of basis used to represent the functional data.

Value

Return glm object plus:

• basis.x: Basis used for fdata or fd covariates.

• basis.b: Basis used for beta parameter estimation.

• beta.l: List of estimated beta parameter of functional covariates.

• data: List that contains the variables in the model.

• formula: Formula.

Note

If the formula only contains a non functional explanatory variables (multivariate covariates), the function compute a standard glm procedure.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Ramsay, James O., and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York.

McCullagh and Nelder (1989), Generalized Linear Models 2nd ed. Chapman and Hall.

Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S, New York: Springer.

See Also

See Also as: predict.fregre.glm and summary.glm.
Alternative method if family=gaussian: fregre.lm.

Examples

## Not run:  
data(tecator)
x=tecator$absorp.fdata
y=tecator$y$Fat
tt=x[["argvals"]]
dataf=as.data.frame(tecator$y)
nbasis.x=11
nbasis.b=7
basis1=create.bspline.basis(rangeval=range(tt),nbasis=nbasis.x)
basis2=create.bspline.basis(rangeval=range(tt),nbasis=nbasis.b)
f=Fat~Protein+x
basis.x=list("x"=basis1)
basis.b=list("x"=basis2)
ldata=list("df"=dataf,"x"=x)
res=fregre.glm(f,family=gaussian(),data=ldata,basis.x=basis.x,
basis.b=basis.b)
summary(res)

## End(Not run)

Variable Selection using Functional Linear Models

Description

Computes functional GLM model between functional covariates (X^1(t_1),\cdots,X^{q}(t_q)) and non functional covariates (Z^1,...,Z^p) with a scalar response Y.

Usage

fregre.glm.vs(
  data = list(),
  y,
  include = "all",
  exclude = "none",
  family = gaussian(),
  weights = NULL,
  basis.x = NULL,
  numbasis.opt = FALSE,
  dcor.min = 0.1,
  alpha = 0.05,
  par.model,
  xydist,
  trace = FALSE
)

Arguments

Details

This function is an extension of the functional generalized spectral additive regression models: fregre.glm where the E[Y|X,Z] is related to the linear prediction \eta via a link function g(\cdot).

E[Y|X,Z]=\eta=g^{-1}(\alpha+\sum_{j=1}^{p}\beta_{j}Z^{j}+\sum_{k=1}^{q}\frac{1}{\sqrt{T_k}}\int_{T_k}{X^{k}(t)\beta_{k}(t)dt})

where Z=\left[ Z^1,\cdots,Z^p \right] are the non functional covariates and X(t)=\left[ X^{1}(t_1),\cdots,X^{q}(t_q) \right] are the functional ones.

Value

Return an object corresponding to the estimated additive mdoel using the selected variables (ame output as thefregre.glm function) and the following elements:

• gof, the goodness of fit for each step of VS algorithm.

• i.predictor, vector with 1 if the variable is selected, 0 otherwise.

• ipredictor, vector with the name of selected variables (in order of selection)

• dcor, the value of distance correlation for each potential covariate and the residual of the model in each step.

Note

If the formula only contains a non functional explanatory variables (multivariate covariates), the function compute a standard glm procedure.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo-de la Fuente manuel.oviedo@udc.es

References

Febrero-Bande, M., Gonz\'alez-Manteiga, W. and Oviedo de la Fuente, M. Variable selection in functional additive regression models, (2018). Computational Statistics, 1-19. DOI: doi:10.1007/s00180-018-0844-5

See Also

See Also as: predict.fregre.glm and summary.glm. Alternative methods: fregre.glm, fregre.glm and fregre.gsam.vs.

Examples

## Not run:  
data(tecator)
x=tecator$absorp.fdata
x1 <- fdata.deriv(x)
x2 <- fdata.deriv(x,nderiv=2)
y=tecator$y$Fat
xcat0 <- cut(rnorm(length(y)),4)
xcat1 <- cut(tecator$y$Protein,4)
xcat2 <- cut(tecator$y$Water,4)
ind <- 1:165
dat <- data.frame("Fat"=y, x1$data, xcat1, xcat2)
ldat <- ldata("df"=dat[ind,],"x"=x[ind,],"x1"=x1[ind,],"x2"=x2[ind,])
# 3 functionals (x,x1,x2), 3 factors (xcat0, xcat1, xcat2)
# and 100 scalars (impact poitns of x1) 

# Time consuming
res.glm0 <- fregre.glm.vs(data=ldat,y="Fat",numbasis.opt=T) # All the covariates
summary(res.glm0)
res.glm0$ipredictors
res.glm0$i.predictor

res.glm1 <- fregre.glm.vs(data=ldat,y="Fat") # All the covariates
summary(res.glm1)
res.glm1$ipredictors
covar <- c("xcat0","xcat1","xcat2","x","x1","x2")
res.glm2 <- fregre.glm.vs(data=ldat, y="Fat", include=covar)
summary(res.glm2)
res.glm2$ipredictors 
res.glm2$i.predictor

res.glm3 <- fregre.glm.vs(data=ldat,y="Fat",
                           basis.x=c("type.basis"="pc","numbasis"=2))
summary(res.glm3)
res.glm3$ipredictors

res.glm4 <- fregre.glm.vs(data=ldat,y="Fat",include=covar,
basis.x=c("type.basis"="pc","numbasis"=5),numbasis.opt=T)
summary(res.glm4)
res.glm4$ipredictors
lpc <- list("x"=create.pc.basis(ldat$x,1:4)
           ,"x1"=create.pc.basis(ldat$x1,1:3)
           ,"x2"=create.pc.basis(ldat$x2,1:4))
res.glm5 <- fregre.glm.vs(data=ldat,y="Fat",basis.x=lpc)
summary(res.glm5)
res.glm5 <- fregre.glm.vs(data=ldat,y="Fat",basis.x=lpc,numbasis.opt=T)
summary(res.glm5)
bsp <- create.fourier.basis(ldat$x$rangeval,7)
lbsp <- list("x"=bsp,"x1"=bsp,"x2"=bsp)
res.glm6 <- fregre.glm.vs(data=ldat,y="Fat",basis.x=lbsp)
summary(res.glm6)
# Prediction like fregre.glm() 
newldat <- ldata("df"=dat[-ind,],"x"=x[-ind,],"x1"=x1[-ind,],
                "x2"=x2[-ind,])
pred.glm1 <- predict(res.glm1,newldat)
pred.glm2 <- predict(res.glm2,newldat)
pred.glm3 <- predict(res.glm3,newldat)
pred.glm4 <- predict(res.glm4,newldat)
pred.glm5 <- predict(res.glm5,newldat)
pred.glm6 <- predict(res.glm6,newldat)
plot(dat[-ind,"Fat"],pred.glm1)
points(dat[-ind,"Fat"],pred.glm2,col=2)
points(dat[-ind,"Fat"],pred.glm3,col=3)
points(dat[-ind,"Fat"],pred.glm4,col=4)
points(dat[-ind,"Fat"],pred.glm5,col=5)
points(dat[-ind,"Fat"],pred.glm6,col=6)
pred2meas(newldat$df$Fat,pred.glm1)
pred2meas(newldat$df$Fat,pred.glm2)
pred2meas(newldat$df$Fat,pred.glm3)
pred2meas(newldat$df$Fat,pred.glm4)
pred2meas(newldat$df$Fat,pred.glm5)
pred2meas(newldat$df$Fat,pred.glm6)

## End(Not run)

Fit Functional Linear Model Using Generalized Least Squares

Description

This function fits a functional linear model using generalized least squares. The errors are allowed to be correlated and/or have unequal variances.

Usage

fregre.gls(
  formula,
  data,
  correlation = NULL,
  basis.x = NULL,
  basis.b = NULL,
  rn,
  lambda,
  weights = NULL,
  subset,
  method = c("REML", "ML"),
  control = list(),
  verbose = FALSE,
  criteria = "GCCV1",
  ...
)

Arguments

Value

An object of class "gls" representing the functional linear model fit. Generic functions such as print, plot, and summary have methods to show the results of the fit.
See glsObject for the components of the fit. The functions resid, coef, and fitted can be used to extract some of its components.
Besides, the class(z) is "gls", "lm", and "fregre.lm" with the following objects:

• sr2: Residual variance.

• Vp: Estimated covariance matrix for the parameters.

• lambda: A roughness penalty.

• basis.x: Basis used for fdata or fd covariates.

• basis.b: Basis used for beta parameter estimation.

• beta.l: List of estimated beta parameter of functional covariates.

• data: List containing the variables in the model.

• formula: Formula used in the adjusted model.

• formula.ini: Formula in call.

• W: Inverse of covariance matrix.

• correlation: See glsObject for the components of the fit.

References

Oviedo de la Fuente, M., Febrero-Bande, M., Pilar Munoz, and Dominguez, A. (2018). Predicting seasonal influenza transmission using functional regression models with temporal dependence. PloS one, 13(4), e0194250. doi:10.1371/journal.pone.0194250

Examples

## Not run:  
data(tecator)
x=tecator$absorp.fdata
x.d2<-fdata.deriv(x,nderiv=)
tt<-x[["argvals"]]
dataf=as.data.frame(tecator$y)

# plot the response
plot(ts(tecator$y$Fat))

nbasis.x=11;nbasis.b=7
basis1=create.bspline.basis(rangeval=range(tt),nbasis=nbasis.x)
basis2=create.bspline.basis(rangeval=range(tt),nbasis=nbasis.b)
basis.x=list("x.d2"=basis1)
basis.b=list("x.d2"=basis2)
ldata=list("df"=dataf,"x.d2"=x.d2)
res.gls=fregre.gls(Fat~x.d2,data=ldata, correlation=corAR1(),
                   basis.x=basis.x,basis.b=basis.b)
summary(res.gls)                   

## End(Not run)

Fitting Functional Generalized Spectral Additive Models

Description

Computes a functional GAM model between a functional covariate (X^1(t_1), \dots, X^{q}(t_q)) and a non-functional covariate (Z^1, ..., Z^p) with a scalar response Y.

This function extends functional generalized linear regression models (fregre.glm) where E[Y|X,Z] is related to the linear predictor \eta via a link function g(\cdot) with integrated smoothness estimation by the smooth functions f(\cdot).

E[Y|X,Z]=\eta=g^{-1}\left(\alpha+\sum_{i=1}^{p}f_{i}(Z^{i})+\sum_{k=1}^{q}\sum_{j=1}^{k_q} f_{j}^{k}(\xi_j^k)\right)

where \xi_j^k is the coefficient of the basis function expansion of X^k; in PCA analysis, \xi_j^k is the score of the j-functional PC of X^k.

Usage

fregre.gsam(
  formula,
  family = gaussian(),
  data = list(),
  weights = NULL,
  basis.x = NULL,
  basis.b = NULL,
  ...
)

Arguments

Details

The smooth functions f(\cdot) can be added to the right-hand side of the formula to specify that the linear predictor depends on smooth functions of predictors using smooth terms s and te as in gam (or linear functionals of these as Z \beta and \langle X(t), \beta \rangle in fregre.glm).

The first item in the data list is called "df" and is a data frame with the response and non-functional explanatory variables, as in gam.

Functional covariates of class fdata or fd are introduced in the following items of the data list.
basis.x is a list of basis functions for representing each functional covariate. The basis object can be created using functions such as create.pc.basis, pca.fd, create.fdata.basis, or create.basis.
basis.b is a list of basis functions for representing each functional beta parameter. If basis.x is a list of functional principal components basis functions (see create.pc.basis or pca.fd), the argument basis.b is ignored.

Value

Return gam object plus:

• basis.x: Basis used for fdata or fd covariates.

• basis.b: Basis used for beta parameter estimation.

• data: List containing the variables in the model.

• formula: Formula used in the model.

• y.pred: Predicted response by cross-validation.

Note

If the formula only contains a non functional explanatory variables (multivariate covariates), the function compute a standard glm procedure.

Author(s)

Manuel Febrero-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Muller HG and Stadtmuller U. (2005). Generalized functional linear models. Ann. Statist.33 774-805.

Wood (2001) mgcv:GAMs and Generalized Ridge Regression for R. R News 1(2):20-25.

Ramsay, James O., and Silverman, Bernard W. (2006), Functional Data Analysis, 2nd ed., Springer, New York.

Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S, New York: Springer.

See Also

See Also as: predict.fregre.gsam and summary.gam.
Alternative methods: fregre.glm and fregre.gkam.

Examples

## Not run: 
data(tecator)
x <- tecator$absorp.fdata
x.d1 <- fdata.deriv(x)
tt <- x[["argvals"]]
dataf <- as.data.frame(tecator$y)
nbasis.x <- 11
nbasis.b <- 5
basis1 <- create.bspline.basis(rangeval=range(tt),nbasis=nbasis.x)
basis2 <- create.bspline.basis(rangeval=range(tt),nbasis=nbasis.b)
f <- Fat ~ s(Protein) + s(x)
basis.x <- list("x"=basis1,"x.d1"=basis1)
basis.b <- list("x"=basis2,"x.d1"=basis2)
ldat <- ldata("df" = dataf, "x" = x , "x.d1" = x.d1)
res <- fregre.gsam(Fat ~ Water + s(Protein) + x + s(x.d1), ldat,
                   family = gaussian(),  basis.x = basis.x,
                   basis.b = basis.b)
summary(res)
pred <- predict(res,ldat)
plot(pred-res$fitted)
pred2 <- predict.gam(res,res$XX)
plot(pred2-res$fitted)
plot(pred2-pred)
res2 <- fregre.gsam(Fat ~ te(Protein, k = 3) + x, data =  ldat,
                     family=gaussian())
summary(res2)

##  dropind basis pc
basis.pc0 <- create.pc.basis(x,c(2,4,7))
basis.pc1 <- create.pc.basis(x.d1,c(1:3))
basis.x <- list("x"=basis.pc0,"x.d1"=basis.pc1)
ldata <- ldata("df"=dataf,"x"=x,"x.d1"=x.d1)  
res.pc <- fregre.gsam(f,data=ldata,family=gaussian(),
          basis.x=basis.x,basis.b=basis.b)
summary(res.pc)
 
##  Binomial family
ldat$df$FatCat <- factor(ifelse(tecator$y$Fat > 20, 1, 0))
res.bin <- fregre.gsam(FatCat ~ Protein + s(x),ldat,family=binomial())
summary(res.bin)
table(ldat$df$FatCat, ifelse(res.bin$fitted.values < 0.5,0,1))

## End(Not run)

Variable Selection using Functional Additive Models

Description

Computes functional GAM model between functional covariates (X^1(t_1),\cdots,X^{q}(t_q)) and non functional covariates (Z^1,...,Z^p) with a scalar response Y.

Usage

fregre.gsam.vs(
  data = list(),
  y,
  include = "all",
  exclude = "none",
  family = gaussian(),
  weights = NULL,
  basis.x = NULL,
  numbasis.opt = FALSE,
  kbs,
  dcor.min = 0.1,
  alpha = 0.05,
  par.model,
  xydist,
  trace = FALSE
)

Arguments

Details

This function is an extension of the functional generalized spectral additive regression models: fregre.gsam where the E[Y|X,Z] is related to the linear prediction \eta via a link function g(\cdot) with integrated smoothness estimation by the smooth functions f(\cdot).

E[Y|X,Z] = \eta = g^{-1}(\alpha + \sum_{i=1}^{p} f_{i}(Z^{i}) + \sum_{k=1}^{q} \sum_{j=1}^{k_q} f_{j}^{k}(\xi_j^k))

where \xi_j^k is the coefficient of the basis function expansion of X^k, (in PCA analysis, \xi_j^k is the score of the j-functional PC of X^k).

The smooth functions f(\cdot) can be added to the right-hand side of the formula to specify that the linear predictor depends on smooth functions of predictors using smooth terms s and te as in gam (or linear functionals of these as Z\beta and \langle X(t), \beta(t) \rangle in fregre.glm).

Value

Return an object corresponding to the estimated additive mdoel using the selected variables (ame output as thefregre.gsam function) and the following elements:

• gof: the goodness of fit for each step of VS algorithm.

• i.predictor: vector with 1 if the variable is selected, 0 otherwise.

• ipredictor: vector with the name of selected variables (in order of selection)

• dcor: the value of distance correlation for each potential covariate and the residual of the model in each step.

Note

If the formula only contains a non functional explanatory variables (multivariate covariates), the function compute a standard gam procedure.

Author(s)

Manuel Feb-Bande, Manuel Oviedo de la Fuente manuel.oviedo@udc.es

References

Febrero-Bande, M., Gonz\'alez-Manteiga, W. and Oviedo de la Fuente, M. Variable selection in functional additive regression models, (2018). Computational Statistics, 1-19. DOI: doi:10.1007/s00180-018-0844-5

See Also

See Also as: predict.fregre.gsam and summary.gam. Alternative methods: fregre.glm, fregre.gsam and fregre.gkam.

Examples

## Not run:  
data(tecator)
x=tecator$absorp.fdata
x1 <- fdata.deriv(x)
x2 <- fdata.deriv(x,nderiv=2)
y=tecator$y$Fat
xcat0 <- cut(rnorm(length(y)),4) 
xcat1 <- cut(tecator$y$Protein,4)
xcat2 <- cut(tecator$y$Water,4)
ind <- 1:165
dat <- data.frame("Fat"=y, x1$data, xcat1, xcat2)
ldat <- ldata("df"=dat[ind,],"x"=x[ind,],"x1"=x1[ind,],"x2"=x2[ind,])
# 3 functionals (x,x1,x2), 3 factors (xcat0, xcat1, xcat2)
# and 100 scalars (impact poitns of x1) 

# Time consuming
res.gam0 <- fregre.gsam.vs(data=ldat,y="Fat"
            ,exclude="x2",numbasis.opt=T) # All the covariates
summary(res.gam0)
res.gam0$ipredictors

res.gam1 <- fregre.gsam.vs(data=ldat,y="Fat") # All the covariates
summary(res.gam1)
res.gam1$ipredictors

covar <- c("xcat0","xcat1","xcat2","x","x1","x2")
res.gam2 <- fregre.gsam.vs(data=ldat, y="Fat", include=covar)
summary(res.gam2)
res.gam2$ipredictors 
res.gam2$i.predictor

res.gam3 <- fregre.gsam.vs(data=ldat,y="Fat",
            basis.x=c("type.basis"="pc","numbasis"=10))
summary(res.gam3)
res.gam3$ipredictors

res.gam4 <- fregre.gsam.vs(data=ldat,y="Fat",include=c("x","x1","x2"),
basis.x=c("type.basis"="pc","numbasis"=5),numbasis.opt=T)
summary(res.gam4)
res.gam4$ipredictors
lpc <- list("x"=create.pc.basis(ldat$x,1:4)
           ,"x1"=create.pc.basis(ldat$x1,1:3)
           ,"x2"=create.pc.basis(ldat$x2,1:12))
res.gam5 <- fregre.gsam.vs(data=ldat,y="Fat",basis.x=lpc)
summary(res.gam5)
res.gam6 <- fregre.gsam.vs(data=ldat,y="Fat",basis.x=lpc,numbasis.opt=T)
summary(res.gam6)
bsp <- create.fourier.basis(ldat$x$rangeval,7)
lbsp <- list("x"=bsp,"x1"=bsp,"x2"=bsp)
res.gam7 <- fregre.gsam.vs(data=ldat,y="Fat",basis.x=lbsp,kbs=4)
summary(res.gam7)
# Prediction like fregre.gsam() 
newldat <- ldata("df"=dat[-ind,],"x"=x[-ind,],"x1"=x1[-ind,],
                "x2"=x2[-ind,])
pred.gam1 <- predict(res.gam1,newldat)
pred.gam2 <- predict(res.gam2,newldat)
pred.gam3 <- predict(res.gam3,newldat)
pred.gam4 <- predict(res.gam4,newldat)
pred.gam5 <- predict(res.gam5,newldat)
pred.gam6 <- predict(res.gam6,newldat)
pred.gam7 <- predict(res.gam7,newldat)
plot(dat[-ind,"Fat"],pred.gam1)
points(dat[-ind,"Fat"],pred.gam2,col=2)
points(dat[-ind,"Fat"],pred.gam3,col=3)
points(dat[-ind,"Fat"],pred.gam4,col=4)
points(dat[-ind,"Fat"],pred.gam5,col=5)
points(dat[-ind,"Fat"],pred.gam6,col=6)
points(dat[-ind,"Fat"],pred.gam7,col=7)
pred2meas(newldat$df$Fat,pred.gam1)
pred2meas(newldat$df$Fat,pred.gam2)
pred2meas(newldat$df$Fat,pred.gam3)
pred2meas(newldat$df$Fat,pred.gam4)
pred2meas(newldat$df$Fat,pred.gam5)
pred2meas(newldat$df$Fat,pred.gam6)
pred2meas(newldat$df$Fat,pred.gam7)

## End(Not run)

Fit of Functional Generalized Least Squares Model Iteratively

Description

This function fits iteratively a functional linear model using generalized least squares. The errors are allowed to be correlated and/or have unequal variances.

1. Begin with a preliminary estimation of \hat{\theta}=\theta_0 (for instance, \theta_0=0). Compute \hat{W}.

2. Estimate b_\Sigma =(Z'\hat{W}Z)^{-1}Z'\hat{W}y

3. Based on the residuals, \hat{e}=\left(y-Zb_\Sigma \right), update \hat{\theta}=\rho\left({\hat{e}}\right) where \rho depends on the dependence structure chosen.

4. Repeats steps 2 and 3 until convergence (small changes in b_\Sigma and/or \hat{\theta}).

Usage

fregre.igls(
  formula,
  data,
  basis.x = NULL,
  basis.b = NULL,
  correlation,
  maxit = 100,
  rn,
  lambda,
  weights = rep(1, n),
  control,
  ...
)

Arguments

Value

An object of class fregre.igls representing the functional linear model fit with temporal dependence errors. Beside, the class(z) is similar to "fregre.lm" plus the following objects:

• corStruct: Fitted AR or ARIMA model.

References

Oviedo de la Fuente, M., Febrero-Bande, M., Pilar Munoz, and Dominguez, A. (2018). Predicting seasonal influenza transmission using functional regression models with temporal dependence. PloS one, 13(4), e0194250. doi:10.1371/journal.pone.0194250

Examples

## Not run:  
data(tecator)
x=tecator$absorp.fdata
x.d2<-fdata.deriv(x,nderiv=)
tt<-x[["argvals"]]
dataf=as.data.frame(tecator$y)
# plot the response
plot(ts(tecator$y$Fat))
ldata=list("df"=dataf,"x.d2"=x.d2)
res.gls=fregre.igls(Fat~x.d2,data=ldata,
correlation=list("cor.ARMA"=list()),
control=list("p"=1)) 
res.gls
res.gls$corStruct

## End(Not run)

Fitting Functional Linear Models

Description

Computes functional regression between functional (and non functional) explanatory variables and scalar response using basis representation.

Usage

fregre.lm(
  formula,
  data,
  basis.x = NULL,
  basis.b = NULL,
  lambda = NULL,
  P = NULL,
  weights = rep(1, n),
  ...
)

Arguments

Details

This section is presented as an extension of the linear regression models: fregre.pc, fregre.pls and fregre.basis. Now, the scalar response Y is estimated by more than one functional covariate X^j(t) and also more than one non functional covariate Z^j. The regression model is given by:

E[Y|X,Z]=\alpha+\sum_{j=1}^{p}\beta_{j}Z^{j}+\sum_{k=1}^{q}\frac{1}{\sqrt{T_k}}\int_{T_k}{X^{k}(t)\beta_{k}(t)dt}

where Z=\left[ Z^1,\cdots,Z^p \right] are the non functional covariates, X(t)=\left[ X^{1}(t_1),\cdots,X^{q}(t_q) \right] are the functional ones and \epsilon are random errors with mean zero , finite variance \sigma^2 and E[X(t)\epsilon]=0.