Datos funcionales

Descripción

Los desarrollos tecnólogicos han hecho posible que los investigadores de muchas áreas dispongan de grandes volúmenes de información para un mismo individuo. Usualmente éstos datos pueden ser representados a través de curvas o en general de funciones.

En el ADF la unidad básica de información es la función completa, más que un conjunto de valores (Ramsay and Dalzell, 1991), es decir, que un dato funcional se puede establecer como la extensión de medidas repetidas cuando la cantidad de mediciones (\(M\)) en un determinado individuo es muy grande.

Ferraty and Vieu (2006) definen una variable aleatoria funcional \(\chi\), como una variable aleatoria que toma valores en un espacio de funciones, es decir, un espacio infinito dimensional (espacio funcional). Una observación \(x\) de la variable aleatoria \(\chi\) se denomina dato funcional.

Un conjunto de datos funcionales \(x_1, x_2, ..., x_n\) es la observación de \(n\) variables funcionales distribuidas como \(\chi\). Un dato funcional \(x_i(t), t ∈ T = [a, b] ⊂ R\), es representado usualmente como un conjunto finito de pares \((t_j, x_{ij})\) con \(t_j\) \(\in\) \(T\), \(j = 1, 2, ..., M\), dónde \(M\) representa la cantidad de puntos en los cuales es observada la variable de interés y \(y_{ij} = \chi_i(t_j)\) (si no existe error observacional) o \(y_{ij} = \chi_i(t_j) + ε_j\) (en caso contrario).

Por ejemplo, consideremos la temperatura promedio mensual de la superficie del océano en grados Celsius, registrada desde Enero de 1982 a Abril de 2022, disponible en:http://www.cpc.noaa.gov/data/indices/sstoi.indices.

Realizamos la lectura de los datos:

library(fda)
## Warning: package 'fda' was built under R version 4.1.3
## Loading required package: splines
## Loading required package: fds
## Warning: package 'fds' was built under R version 4.1.3
## Loading required package: rainbow
## Warning: package 'rainbow' was built under R version 4.1.3
## Loading required package: MASS
## Loading required package: pcaPP
## Warning: package 'pcaPP' was built under R version 4.1.3
## Loading required package: RCurl
## Warning: package 'RCurl' was built under R version 4.1.2
## Loading required package: deSolve
## Warning: package 'deSolve' was built under R version 4.1.3
## 
## Attaching package: 'fda'
## The following object is masked from 'package:graphics':
## 
##     matplot
library(fda.usc)
## Warning: package 'fda.usc' was built under R version 4.1.3
## Loading required package: mgcv
## Warning: package 'mgcv' was built under R version 4.1.1
## Loading required package: nlme
## This is mgcv 1.8-37. For overview type 'help("mgcv-package")'.
## ----------------------------------------------------------------------------------
##  Functional Data Analysis and Utilities for Statistical Computing
##  fda.usc version 2.0.2 (built on 2020-02-17) is now loaded
##  fda.usc is running sequentially usign foreach package
##  Please, execute ops.fda.usc() once to run in local parallel mode
##  Deprecated functions: min.basis, min.np, anova.hetero, anova.onefactor, anova.RPm
##  New functions: optim.basis, optim.np, fanova.hetero, fanova.onefactor, fanova.RPm
## ----------------------------------------------------------------------------------
library(rainbow)
library(MASS)
library(xtable)
setwd("C:/Users/jpari/Dropbox/FDA")
Datos =read.delim2("sstoi.indices1.txt",header=T,dec=".", sep = "\t")
summary(Datos)
##        YR            MON            NINO1.2           ANOM         
##  Min.   :1982   Min.   : 1.000   Min.   :19.06   Min.   :-1.90000  
##  1st Qu.:1992   1st Qu.: 3.000   1st Qu.:21.23   1st Qu.:-0.73250  
##  Median :2002   Median : 6.000   Median :23.14   Median :-0.23500  
##  Mean   :2002   Mean   : 6.467   Mean   :23.25   Mean   :-0.05118  
##  3rd Qu.:2012   3rd Qu.: 9.000   3rd Qu.:25.22   3rd Qu.: 0.44000  
##  Max.   :2022   Max.   :12.000   Max.   :28.51   Max.   : 4.03000  
##      NINO3           ANOM.1             NINO4           ANOM.2        
##  Min.   :23.38   Min.   :-2.16000   Min.   :26.36   Min.   :-1.87000  
##  1st Qu.:25.00   1st Qu.:-0.65000   1st Qu.:28.02   1st Qu.:-0.55250  
##  Median :25.93   Median :-0.15000   Median :28.57   Median :-0.00500  
##  Mean   :25.97   Mean   :-0.06033   Mean   :28.46   Mean   :-0.08946  
##  3rd Qu.:26.89   3rd Qu.: 0.42000   3rd Qu.:28.98   3rd Qu.: 0.38000  
##  Max.   :28.81   Max.   : 3.07000   Max.   :30.22   Max.   : 1.55000  
##     NINO3.4          ANOM.3        
##  Min.   :24.56   Min.   :-2.22000  
##  1st Qu.:26.35   1st Qu.:-0.65250  
##  Median :27.07   Median :-0.11000  
##  Mean   :27.02   Mean   :-0.06878  
##  3rd Qu.:27.69   3rd Qu.: 0.45000  
##  Max.   :29.54   Max.   : 2.72000

Reorganizamos la información para que cada año sea una curva. No tenemos en cuenta el año 2022 por no disponer de la información de todo el año.

temp <- matrix(Datos$NINO1.2[-c(481:484)], nrow=12, 40)
colnames(temp)=c(1982:2021)

Construimos las gráficas para cada año

month<-1:12
plot(month, temp[,1], type="b", ylim=c(18, 30), xlab="Mes", ylab="temperatura")
for(i in 2:ncol(temp))
lines(month, temp[,i], type="b", col=i)

Suavizamiento de curvas

Una herramienta no paramétrica de mucha utilidad en el ADF es el suavizado de curvas a través de funciones. El procedimiento consiste en aproximar las funciones del espacio considerado a través de (Ramsay and Silverman, 2005):

\[\begin{equation} \chi(t) = \sum_{l=1}^K c_lB_l(t) = c'B(t) \end{equation}\] dónde K es el número de funciones de la base.

En la literatura exisste varias funciones base que permiten desarrollar la expansión. La selección del tipo de funciones de la base depende de las caracteríticas que cumpla que fenómeno de estudio. Por ejemplo, las series de Fourier se utilizan para funciones con comportamientos cíclicos o periódicos y las bases monomiales se utilizan cuando existen tendencias simples que pueden ser ajustadas mediante líneas rectas, polinomios cuadráticos, polinomios de orden superior, etc. Las bases más usuales son aquellas basadas en splines debido a que son más flexibles y se ajustan de mejor manera a diferentes comportamientos. Dentro de estas se encuentran B-splines, M-splines, I-splines, y fucniones de potencia truncadas.

Bases monomiales

Las bases monomiales requieren el dominio y el número de base. Por ejemplo, la base monomial con K=6 funciones de base definidas en el intervalo [0,1] se puede construir con:

library(fda)
bbasis_obj = create.monomial.basis(rangeval=c(0,1), nbasis = 6)

Esto devolverá una salida de “funciones”. Para evaluar las bases en una cuadrícula de puntos \(s\), debemos:

library(fda)
x <- seq(0,1,length.out=100)
bbasisevals <- eval.basis(x, bbasis_obj)
# dim(basisevals)
matplot(x, bbasisevals, type='l', lty=1, col=rainbow(6),
        xlab="x", ylab="bases", 
        main="base monomial con k = 6")

Bases de Fourier

Para utilizarlas bases de Fourier se requiere que se defina el dominio, el período de oscilación y el número de funciones de base.

fbasis_obj <- create.fourier.basis(rangeval=c(0,1), 
                                   nbasis=65, period = 1)
fbasisevals <- eval.basis(x, fbasis_obj)
matplot(x, fbasisevals[, 1:3], type='l', lty=1, col=rainbow(3),
        xlab="x", ylab="bases", 
        main="Primeras tres bases de Fourier")

B-splines

Las bases B-spline requieren el dominio, el número de funciones de la base y el orden.

bsbasis_obj <- create.bspline.basis(rangeval=c(0,1),
                                    nbasis=10, norder=4)
bsbasisevals <- eval.basis(x, bsbasis_obj)
matplot(x, bsbasisevals, type='l', lty=1, col=rainbow(15),
        xlab="x", ylab="bases", 
        main="B-spline cubica con K = 10")

Otras bases

El paquete también se puede usar para construir otros tipos de bases. Para revisar la lista de bases disponibles se puede utilizar $ ?create. + tab$

Suavizamiento de la temperatura

La temperatura es una variable que se puede tomar en intervalos más pequeños (diaria, cada hora, etc). Es decir, al realizar la gráfica anterior estamos suponiendo que podemos discretizar la variable de temperatura a través de la medición mensual.La primera tarea es convertir estos valores discretos en una función que puede tomar valores en cualquier valor de argumento deseado \(x_i(t)\). Si se supone que estas observaciones no presentan un término de error asociado, el proceso de generación de la función (curva) será un proceso de interpolación, pero si tienen algún error observacional que necesita modelarse, entonces la conversión de datos (finitos) a funciones (que teóricamente se pueden evaluar en un número infinito de puntos) se enmarca en los procesos de suavizamiento.

En consecuencia, requerimos una estrategia para construir funciones con parámetros que sean fáciles de estimar y que puedan ajustarse casi a cualquier característica de la curva. Por otro lado, no queremos usar más parámetros de los que necesitamos, ya que hacerlo aumentaría en gran medida el tiempo de cálculo y complicar los análisis. Para el caso de la temperatura se utiliza una B-spline de orden 4.

BSpl <- create.bspline.basis(norder=4, breaks=seq(0,12,length=5))
plot(BSpl)

Ftemp <- Data2fd(temp, basisobj=BSpl)
## 'y' is missing, using 'argvals'
## 'argvals' is missing;  using seq( 0 ,  12 , length= 12 )
plot(Ftemp, main="Datos suavizados de la Temperatura", xlab = "Mes", ylab = "Temperatura oC",
ylim=c(18,30))

## [1] "done"

En el proceso de suavizamiento es necesario tener en cuenta la rugosidad como un aspecto fundamental. El método utilizado para aproximar funciones mediante una base de funciones B-splines se basa en minimizar el sistema:

\[\begin{equation} GCV = \min_c \sum_{j =1}^M (y_j - S(t_j))^2 + \lambda\int_t (S''(t)dt) \end{equation}\]

con \(\lambda\) un parámetro de penalización para disminuir la variabilidad del ajuste. El número de funciones base (K) y el coeficiente \(\lambda\) son estimados a través de procedimientos de validación cruzada (Ramsay and Silverman, 2005). Con el siguiente programa se puede establecer la estimación del parámetro de suavizamiento óptimo.

loglam = seq(0,0.05,0.001)
nlam = length(loglam)
dfsave = rep(NA,nlam)
gcvsave = rep(NA,nlam)
for (ilam in 1:nlam) {
lambda = loglam[ilam]
fdParobj = fdPar(BSpl, Lfdobj=NULL, lambda= lambda)
smoothlist = smooth.basis(1:12, temp, fdParobj)
dfsave[ilam] = smoothlist$df
gcvsave[ilam] = sum(smoothlist$gcv)
}
plot(loglam,gcvsave,xlab=expression(lambda),ylab=expression(GCV(lambda)),
main="Parámetros de suavizamiento versus GCV",type="b",cex=0.7)

Tomando como referencia los 40 añoos reportados (de 1982 a 2021) en la base de datos de interés, se evidencia que el punto que minimiza la estadística GCV se encuentra alrededor de \(\nu\) = 0.017 con valores muy superiores antes de 0.01 y después de 0.04. Se destaca que dicho parámetro minimiza la suma del criterio GCV de todas las 40 curvas que se encuentran dentro de la muestra, es decir, el valor de dicho parámetro resulta ser el óptimo en relación con el proceso de suavizamiento de todas las curvas dentro de la muestra.

Estadística descriptiva

Dado un conjunto de datos funcionales \(x_1, x_2, ..., x_n\) definidos en \(t ⊂ T ∈ R\), las correspondientes funciones descriptivas están dadas por las expresiones (Ramsay and Silverman, 2005):

\(Media:\bar{\chi}(t) = \frac{\sum_{i=1}^n \chi_i(t)}{n}\)

\(Varianza: Var(\chi(t)) = \frac{\sum_{i=1}^n (\chi_i(t) - \bar{\chi}(t))^2}{n}\)

\(Desviación estándar: \sigma(\chi(t)) = \sqrt{var(\chi(t)}\)

Así, se puede concluir que las estadísticas descriptivas univariadas clásicas se aplican igualmente cuando se tienen datos funcionales. Sin embargo se resalta que en este caso, los objetos calculados corresponden a curvas.

meanfdh <- mean.fd(Ftemp)
varfdh <-var.fd(Ftemp)
stdvfdh <- stddev.fd(Ftemp)
plot(Ftemp,col=8, lty=1, xlab = "Mes", ylab = "Temperatura oC")
## [1] "done"
lines(meanfdh,col=2,lwd=2)

par(mfrow=c(1,2))
plot(varfdh, main=" Superficie de varianza", xlab = "t", ylab = "s")

plot(stdvfdh, main="Desviación estándar", xlab = "Mes", ylab = "Temperatura oC")

## [1] "done"