Código
library(R.matlab)
library(fda)
library(plotly)
library(dplyr)
library(depthTools)
library(mrfDepth)
library(ggplot2)
library(dplyr)
library(tidyr)
library(robustbase)
library(Matrix)
library(kableExtra)
library(fdaoutlier)Taller 1 - Análisis de datos funcionales.
1 Descripción de datos
Fluorescencia
La fluorescencia es la propiedad de algunos átomos y moléculas de absorber luz a una longitud de onda determinada (excitación) seguido por la emisión de luz de corta duración a una longitud de onda más larga.
La fluorescencia implica una fuente de luz externa para excitar la muestra a una longitud de onda en particular. Cuando se excita a la longitud de onda adecuada, la molécula pasa de un estado básico a otro excitado. A medida que la molécula vuelve al estado fundamental, se libera energía en forma de calor (pérdida de energía) y luz a una longitud de onda diferente de menor energía.
Recolección de los datos
Se recolectaron muestras de azúcar durante tres meses de operación en una planta azucarera ubicada en Escandinavia. Las muestras se tomaron cada ocho horas, lo que corresponde a tres mediciones diarias, para un total de 268 muestras a lo largo del periodo de estudio.
El azúcar se disolvió en agua en una proporción de 2.25 g por 15 ml, y la solución resultante se analizó utilizando un espectrofluorómetro PE LS50B.
Cada muestra fue excitada a siete longitudes de onda (230, 240, 255, 290, 305, 325 y 340 nm). Para cada una de estas excitaciones, el equipo registró la intensidad de emisión en 571 longitudes de onda. Los espectros de emisión se midieron en el rango de 275 a 560 nm, con intervalos de 0,5 nm.
Almacenamiento de los datos
Los datos de fluorescencia se almacenan en una matriz de 268 x 3997. Las primeras 571 columnas corresponden a las mediciones de emisión para la primera longitud de onda de excitación (340 nm). Las últimas 571 columnas corresponden a las mediciones de emisión para la última longitud de onda de excitación (230 nm)
Asi se tienen las curvas espectrométricas asociadas a la producción de azúcar corresponden a 7 longitudes de onda de excitación \(Y_1 = 230\,\text{nm},\; Y_2 = 240\,\text{nm},\; Y_3 = 255\,\text{nm},\; Y_4 = 290 \text{nm},\; Y_5 = 305\,\text{nm},\)
\(Y_6 = 325\,\text{nm},\; Y_7 = 340\,\text{nm}\).
Las observaciones correspondientes a cada longitud de excitación representan un proceso funcional continuo \(Y_i\) , así que el proceso de producción de azúcar se puede considerar como un proceso funcional multivariado \(\mathbf{Y}\) está conformado por \(\mathbf{Y}=(Y_1, Y_2, ..., Y_7)\).
Código
sugar <- sugar <- readMat("C:\\Users\\sofia\\OneDrive\\Documentos\\2026-1\\FDA\\data.mat")
sugar_fluo <- sugar$X
dim(sugar_fluo) [1] 268 3997Código
t <-sugar$EmAx
wvl <- array(NA, dim = c(571, 268, 7))
for (k in 1:7) {
start <- 571 * (k- 1) + 1
end <- 571 * k
wvl[, , 8- k] <- t(sugar_fluo[,start:end])
}
nms <- c("340nm", "325nm", "305nm", "290nm", "255nm", "240nm", "230nm")
nms <- rev(nms)
dimnames(wvl) <- list(NULL, NULL, nms)2 Visualización de los espectros de emisión
3 Suavizado de las curva
El objetivo es ajustar una función a un conjunto de observaciones discretizadas \(y_j\), con \(j = 1, \ldots, n\), bajo \[ y_j = x(t_j) + \varepsilon_j \]
Donde \(x\) es una función desconocida y \(\varepsilon_j\) es el término de error.
La función \(x\) para cada tiempo, se aproxima mediante:
\[ \hat{x}(t) \approx \sum_{k=1}^{K} c_k \phi_k(t) \]
Las funciones \(\phi_k\) constituyen la base funcional utilizada en la aproximación y su forma específica depende del tipo de base seleccionada (por ejemplo, Fourier o B-splines).
Los coeficientes \(c_k\) representan los pesos asociados a cada función base dentro de la combinación lineal y se estiman mediante mínimos cuadrados o mínimos cuadrados penalizados. Finalmente, \(K\) es el número total de funciones base empleadas en la aproximación y determina, en parte, la flexibilidad del modelo.
Suavizado splines
La metodología de suavizamiento mediante splines también se basa en la representación de la función como una combinación lineal de funciones base. Sin embargo, en este caso, los coeficientes de dicha combinación se estiman minimizando una suma de cuadrados del error a la que se le incorpora un término de penalización asociado a la curvatura de la función, controlado por el parámetro \(\lambda\).
\(\mathrm{PENSSE}_m(y \mid c) = \mathrm{SSE}(y \mid c) + \lambda \mathrm{PEN}_m(x)\)
Este parámetro regula el equilibrio entre ajuste y suavidad. Generalmente, \(\lambda\) se selecciona mediante validación cruzada generalizada (GCV).
En esta metodología, a diferencia de la regresión splines, el número de funciones base no es el elemento más determinante; el papel central lo desempeña el valor de \(\lambda\), ya que es el que define principalmente el grado de suavidad de la función.
Justificación del suavizado
Para el suavizado de las curvas se empleará suaviado splines, utilizando bases B-splines. Esta elección se justifica a partir de la visualización de los datos, donde se puede observar en Figura 1 que los espectros de emisión no presentan un comportamiento periódico. En este contexto, los B-splines son una base adecuada.
En cuanto al número de funciones base, se realizaron pruebas con distintos valores de \(K\) (por ejemplo, valores cercanos a 70). No obstante, se observó que algunos de estos ajustes producían valores negativos en los espectros de emisión, lo cual no es consistente con la naturaleza de los datos. Por esta razón, se utilizo un número de funciones base acorde con la cantidad total de puntos observados (571) , para cada uno de los siete procesos analizados.
Finalmente, la selección del parámetro de suavizado \(\lambda\) se llevó a cabo mediante la exploración de una grilla de valores candidatos, con el objetivo de identificar aquel que minimiza el criterio de validación cruzada generalizada (GCV). Esta estrategia de selección del parámetro de suavizado sigue lo propuesto en Ramsay, Hooker, y Graves (2009).
Código
x <- seq(275, 560, by = 0.5)
onda_340 <- as.matrix(t(sugar_fluo[1:268,1:571]))
onda_325 <- as.matrix(t(sugar_fluo[1:268,572:1142]))
onda_305 <- as.matrix(t(sugar_fluo[1:268,1143:1713]))
onda_290 <- as.matrix(t(sugar_fluo[1:268,1714:2284]))
onda_255 <- as.matrix(t(sugar_fluo[1:268,2285:2855]))
onda_240 <- as.matrix(t(sugar_fluo[1:268,2856:3426]))
onda_230 <- as.matrix(t(sugar_fluo[1:268,3427:3997]))Código
ondas_lista <- list(
"340" = onda_340, "325" = onda_325, "305" = onda_305,
"290" = onda_290, "255" = onda_255, "240" = onda_240,
"230" = onda_230
)
resultados_suavizados <- list()
lambdas <- 10^seq(-3, 4, 0.05)
BSpl <- create.bspline.basis(
rangeval = c(min(x), max(x)),
nbasis = 571
)
for (nombre in names(ondas_lista)) {
data_actual <- ondas_lista[[nombre]]
gcv_values <- numeric(length(lambdas))
for(i in 1:length(lambdas)){
fdParobj <- fdPar(BSpl, Lfdobj=2, lambda=lambdas[i])
Fsugar_tmp <- smooth.basis(x, data_actual, fdParobj)
gcv_values[i] <- sum(Fsugar_tmp$gcv)
}
mejor_lambda <- lambdas[which.min(gcv_values)]
final_fdPar <- fdPar(BSpl, Lfdobj=2, lambda=mejor_lambda)
resultados_suavizados[[nombre]] <- smooth.basis(x, data_actual, final_fdPar)
cat("Mejor lambda:", mejor_lambda, "\n")
}
saveRDS(resultados_suavizados, file = "suavizados.rds")| Parámetro | Valor Seleccionado |
|---|---|
| \(\lambda_1\) | 0.28183829 |
| \(\lambda_2\) | 0.25118864 |
| \(\lambda_3\) | 0.28183829 |
| \(\lambda_4\) | 0.06309573 |
| \(\lambda_5\) | 0.06309573 |
| \(\lambda_6\) | 0.1 |
| \(\lambda_7\) | 0.05623413 |
Tabla 1: Valores seleccionados del parámetro de suavizado mediante GCV.
Código
setwd("C:\\Users\\sofia\\OneDrive\\Documentos\\2026-1\\FDA")
resultados_suavizados <- readRDS("suavizados.rds")4 Usando el data set, seleccione un proceso y encuentre:
Se considerará para los puntos a continuación el proceso \(Y_2\) que corresponde a las curvas espectrométricas asociadas a la producción de azúcar de la segunda longitud de onda de excitación (240nm).
4.1 Función media
Sea \(x_i\), \(i=1, \ldots, n\), una muestra de curvas o funciones ajustadas a los datos. La función media muestral en \(t\), se definen como:
\[ \bar{x}(t) = \frac{\sum_{i=1}^n x_i(t)}{n} \]
Usando la función mean.fd del paquete fda se presenta a continuación la función media muestral.
Código
fd_240 <- resultados_suavizados$`240`$fd
media_240 <- mean.fd(fd_240)
plot(fd_240,
col = "#CCCCFF",
lty = 1,
xlab = "Longitud de onda (nm)",
ylab = "Intensidad")[1] "done"Código
lines(media_240, col = "#000033", lwd = 3)
legend("topright",
legend = "Función media",
col = "#000033",
lwd = c(1, 3))
En Figura 2 se presentan las curvas espectrométricas asociadas a la segunda longitud de onda de excitación. Las líneas azules representan las curvas individuales, mientras que la línea azul oscuro corresponde a la función media muestral. Se observa un aumento de la intensidad hasta aproximadamente 350 nm, seguido de un comportamiento casi constante hasta alrededor de 450 nm. A partir de este punto, la intensidad comienza a descender, lo que indica que para longitudes de onda menores y mayores la intensidad es menor.
4.2 Función media recortada al 10%
La media recortada al \(100\alpha \%\) se calcula de manera análoga a la media funcional, con la diferencia de que se excluyen las \(100\alpha \%\) funciones menos centrales.
Functional band depth
Siguiendo el enfoque propuesto por López-Pintado y Romo (2009). Sea \(x_1, \dots, x_n\) una colección de funciones reales, en el espacio \(C(I)\) de funciones reales continuas sobre el intervalo compacto \(I\). El gráfico de una función \(x\) es el subconjunto del plano \(G(x) = \{(t, x(t)) : t \in I\}\). La banda en \(\mathbb{R}^2\) delimitada por las curvas \(x_{i_1}, \dots, x_{i_k}\) es:
\[ \begin{equation}B(x_{i_1}, x_{i_2}, \dots, x_{i_k}) = \left\{ (t, y) : t \in I, \min_{r=1, \dots, k} x_{i_r}(t) \leq y \leq \max_{r=1, \dots, k} x_{i_r}(t) \right\}\end{equation} \]
Para cualquier función \(x\) en \(x_1, \dots, x_n\), la cantidad
\[ \begin{equation}S_n^{(j)}(x) = \binom{n}{j}^{-1} \sum_{1 \leq i_1 < i_2 < \dots < i_j \leq n} I\{G(x) \subset B(x_{i_1}, x_{i_2}, \dots, x_{i_j})\}, \quad j \geq 2,\end{equation} \]
expresa la proporción de bandas \(B(x_{i_1}, x_{i_2}, \dots, x_{i_j})\) determinadas por \(j\) curvas diferentes \(x_{i_1}, x_{i_2}, \dots, x_{i_j}\) que contienen el gráfico de \(x\).
Para las funciones \(x_1, \dots, x_n\), la profundidad de banda (BD) de cualquiera de estas curvas \(x\) es
\[ \begin{equation}S_{n,J}(x) = \sum_{j=2}^{J} S_n^{(j)}(x), \quad J \geq 2.\end{equation} \]
De esta manera, es posible ordenar las curvas según su nivel de centralidad, de modo que aquellas con mayores valores de profundidad de banda (BD) se consideran las más centrales dentro del conjunto de datos.
Tambien se puede introducir una definición más flexible donde \(GS_n^{(j)}(x)\) es la proporción de coordenadas de \(x\) dentro del intervalo establecido por \(j\) puntos diferentes de la muestra:
\[ \begin{equation}GS_n^{(j)}(x) =\binom{n}{j}^{-1}\sum_{1 \leq i_1 < \dots < i_j \leq n}\frac{1}{d}\sum_{k=1}^{d}I \left\{\min \{ x_{i_1}(k), \dots, x_{i_j}(k) \}\leq x(k) \leq\max \{ x_{i_1}(k), \dots, x_{i_j}(k) \}\right\}.\end{equation} \]
Definiendo este orden funcional, sean \(x_{[1]}, \dots, x_{[n]}\) los estadísticos de orden, siendo \(x_{[1]}\) la curva más profunda (la más central) y \(x_{[n]}\) la curva más alejada.
La media recortada al \(100\alpha\%\), se calcula
\[ \overline{x}_{\alpha}(t) = \sum_{i=1}^{n - \lfloor n\alpha \rfloor} \frac{x_{[i]}(t)}{n - \lfloor n\alpha \rfloor} \]
De esta forma, si se quiere encontrar la media recortada al 10%, entonces con \(\alpha = 0.1\) y \(N=268\), \(N - \lfloor N\alpha \rfloor = 242\).
Usando la función MBD del paquete depthTools se calcula esta última propuesta de profundidad de banda. Posteriormente, se presenta la visualización de las curvas junto con la media recortada.
Código
n = 267
alpha <- 0.1
upto <- floor(n * alpha)
mbd <- MBD(t(onda_240), plotting = FALSE)
depths <- data.frame(idx = seq_along(mbd$MBD),
MBD = as.numeric(mbd$MBD))
depths <- depths[order(depths$MBD, decreasing = F), ]
ind <- depths$idx[1:upto]
fd_240_10 <- fd_240
fd_240_10$coefs <- fd_240_10$coefs[, -c(ind)]
media_240_10 <- mean.fd(fd_240_10)
media_240 <- mean.fd(fd_240)
plot(fd_240,
col = "#CCCCFF",
lty = 1,
xlab = "Longitud de onda (nm)",
ylab = "Intensidad")[1] "done"Código
lines(media_240, col = "#000033", lwd = 3)
lines(media_240_10, col = "#D9D9D9", lwd = 3)
legend("topright",
legend = c("Función media","Función media recortada"),
col = c("#000033","#D9D9D9"),
lwd = c(1, 3))
4.3 Función de varianza
Sea \(x_1, ..., x_n\) una muestra de curvas sobre un intervalo \([a,b]\). Luego, la función de varianza en \(t\) está dada por:
\[s(t) = \frac{1}{n-1}\sum_{i=1}^n (x_i (t) - \bar{x}(t))^2\]
Esta función se ilustra en Figura 4.
Código
varianza_funcional <- var.fd(fd_240)
cov_mat2 <- eval.bifd(x, x, varianza_funcional)
var_diag <- diag(cov_mat2)
plot(x, var_diag, type="l", lwd=1,
xlab="Longitud de Onda de Excitación (nm)", ylab="Varianza", col="red", main = "Función Varianza") 
4.4 Función de covarianza
La función de covarianza bivariada \(\sigma(s,t)\) especifica la covarianza entre los valores de las curvas \(x_i(s)\) y \(x_i(t)\) en los tiempos \(s\) y \(t\), respectivamente. Esta se estima mediante
\(v(s,t) = (n-1)^{-1} \sum_{i=1}^{n} [x_i(s) - \bar{x}(s)][x_i(t) - \bar{x}(t)]\)
Código
sugarvar.bifd = var.fd(resultados_suavizados$`240`$fd)
sugarvar_mat = eval.bifd(x, x, sugarvar.bifd)
escala_colores <- list(
list(0, "#000033"),
list(0.02, "#0000FF"),
list(0.05, "#FFFFFF"),
list(0.1, "#FFFF00"),
list(0.3, "#FF7F00"),
list(0.6, "#FF0000"),
list(1, "#4B0082")
)
p <- plot_ly(x = ~x, y = ~x, z = ~sugarvar_mat) %>%
add_surface(
colorscale = escala_colores,
colorbar = list(title = "Covarianza/"),
contours = list(
z = list(
show = TRUE,
usecolormap = TRUE,
highlightcolor = "#ff0000",
project = list(z = TRUE)
)
)
) %>%
layout(
title = list(text = "Superficie de Covarianza (Onda 240)"),
scene = list(
xaxis = list(title = "Longitud de onda (nm)"),
yaxis = list(title = "Longitud de onda (nm)"),
zaxis = list(
title = "Valor",
range = c(0,334),
tickformat = ".0f"
)
)
)
pEn Figura 5 se observa que para longitudes de onda pequeñas y grandes los valores son cercanos a cero, lo que indica una menor variabilidad conjunta de las intensidades en esas regiones del espectro. Esto es consistente con el comportamiento observado en las curvas, donde las funciones tienden a concentrarse alrededor de la media en los extremos del dominio.
En contraste, en la región central del espectro, aproximadamente entre 350 y 500 nm, la covarianza aumenta considerablemente y alcanza su valor máximo cerca de 400 nm. Esto sugiere que las intensidades en esas longitudes de onda presentan una mayor variabilidad entre las muestras y además tienden a variar conjuntamente.
En consecuencia, esta región concentra la mayor parte de la variabilidad del proceso de fluorescencia, lo cual también se refleja en la mayor dispersión de las curvas alrededor de la función media.
4.5 Función mediana
La funcion mediana muestral es la curva con la mayor profundidad. Es decir, \[ \hat{m}_{n,J} = \underset{x \in \{x_1, \ldots, x_n\}}{\arg \min} MBD_n(x) \]
Esta función se ilustra en la Figura 6.
Código
depths <- depths[order(depths$MBD, decreasing = T), ]
mediana <- depths$idx[1]
fd_240_med <- fd_240
fd_240_med$coefs <- fd_240_med$coefs[, mediana, drop = FALSE]
media_240 <- mean.fd(fd_240)
plot(fd_240,
col = "#CCCCFF",
lty = 1,
xlab = "Longitud de onda (nm)",
ylab = "Intensidad")[1] "done"Código
lines(x, eval.fd(x, media_240), col="#000033", lwd=3)
lines(x, eval.fd(x,fd_240_med), col = "#CC0000", lwd = 3)
legend("topright",
legend = c("Media", "Mediana funcional"),
col = c("#000033", "#CC0000"),
lwd = c(1,3))
4.6 Funciones quantiles 90 y 95
Para hallar las funciones se utiliza la relación \(p=\mid 2q-1 \mid\), donde \(p\) es el porcentaje de curvas más profundas y \(q\) representa el cuantil.
4.6.1 Función cuantil 90
En este caso, se sabe que \(q=0.9\), luego \(p=0.8\). Así, la función cuantil 90 se ilustra en la Figura 7.
Código
curvas_eval <- eval.fd(x, fd_240)
mbd_values <- modified_band_depth(dt = t(curvas_eval))
# Con la relación p=|2q-1|
p = abs(2*0.9-1)
n = 268
k <- ceiling(p * n) # cantidad que corresponde al 80% de curvas más profundas
orden <- order(mbd_values, decreasing = TRUE)
indices80 <- orden[1:k]
curvas80 <- fd_240[indices80]
class(curvas80)[1] "fd"Código
tgrid <- seq(min(x), max(x), length = 200)
Y80 <- eval.fd(x, curvas80)
curva_min <- apply(Y80, 1, min)
# Función cuantil 90
curva_max <- apply(Y80, 1, max)
# Graficando todas las curvas
plot(fd_240, col = "lightgray", xlab="Longitud de Onda de Excitación (nm)", ylab="Intensidad de Fluorescencia")[1] "done"Código
# Graficando q1 = 0.9
lines(x, curva_max, col="tomato", lwd=2)
legend("topright",
legend = bquote("Función Cuantil 90"),
col = "tomato",
lwd = 3,
bty = "n",
cex = 0.9)
4.6.2 Función cuantil 95
En este caso, se sabe que \(q=0.95\), luego \(p=0.9\). Así, la función cuantil 90 se ilustra en la Figura 8.
Código
# Con la relación p=|2q-1|
# Con la relación p=|2q-1|
p = abs(2*0.95-1)
k <- ceiling(p * n) # cantidad que corresponde al 90%
indices90 <- orden[1:k]
curvas90 <- fd_240[indices90]
class(curvas90)[1] "fd"Código
Y90 <- eval.fd(x, curvas90)
curva_min <- apply(Y90, 1, min)
# Función cuantil 95
curva_max <- apply(Y90, 1, max)
# Graficando todas las curvas
plot(fd_240, col = "lightgray", xlab="Longitud de Onda de Excitación (nm)", ylab="Intensidad de Fluorescencia")[1] "done"Código
# Graficando q1 = 0.95 y q2 = 0.1
lines(x, curva_max, col="tomato", lwd=2)
legend("topright",
legend = bquote("Función Cuantil 95"),
col = "tomato",
lwd = 3,
bty = "n",
cex = 0.9)
4.7 Región central del 0.75
Sea \(y_{[i]}\) la curva de la muestra asociada con el \(i\)-ésimo valor más grande de profundidad de banda. Considerando a \(y_{[1]}, \dots, y_{[n]}\) como estadísticos de orden, siendo \(y_{[1]}\) la curva más profunda (la más central) y \(y_{[n]}\) la curva más alejada.
La región central del \(0.75\%\) de la muestra se expresa por
\[C_{0.75} = \left\{ (t,y(t)): \min_{r=1,\dots, \lceil 0.75n \rceil } y_{[r]}(t) \leq y(t) \leq \max_{r=1, \dots, \lceil 0.75n \rceil } y_{[r]}(t) \right\}\]
Partiendo de la relación \(p=|2q-1|\)
\[ \begin{equation} \begin{aligned} p &= 0.75 \end{aligned} \end{equation} \] \[ \begin{equation} \begin{aligned} 2q - 1 &= 0.75 \\ 2q &= 1.75 \\ q &= \frac{1.75}{2} \\ q &= 0.875 \end{aligned} \end{equation} \] \[ \begin{equation} \begin{aligned} -2q + 1 &= 0.75 \\ 2q &= 1 - 0.75 \\ q &= \frac{0.25}{2} \\ q &= 0.125 \end{aligned} \end{equation} \]
Usando la función MBD del paquete depthTools
Código
res_mbd <- MBD(t(onda_240), plot = FALSE)
valores_mbd <- res_mbd$MBD
orden <- order(valores_mbd, decreasing = TRUE)
k <- ceiling(0.75 * length(valores_mbd))
central <- orden[1:k]
curvas_centrales <- t(onda_240)[central, ]
lim_inf <- apply(curvas_centrales, 2, min)
lim_sup <- apply(curvas_centrales, 2, max)
curvas <- eval.fd(x, resultados_suavizados$`240`$fd)
matplot(x, curvas, type="l", col="#CCCCFF", lty=1,
xlab="Longitud de onda (nm)",
ylab="Intensidad")
polygon(c(x, rev(x)),
c(lim_inf, rev(lim_sup)),
col=rgb(0,0,1,0.2), border=NA)
lines(x, lim_inf, lwd=2, col="#000033")
lines(x, lim_sup, lwd=2, col="#000033")
En la Figura 9 se observa que la región sombreada en azul representa la región central del \(75\%\), construida a partir de las curvas con mayor profundidad de banda. Las líneas negras gruesas que delimitan esta zona corresponden a la envolvente superior e inferior, definida matemáticamente como el valor mínimo y el valor máximo en cada punto \(t\) entre las \(0.75n\) curvas más profundas, excluyendo así el \(25\%\) de las curvas más periféricas
4.8 Outliers usando el boxplot functional
Para poder construir el boxplot funcional, primero se debe hallar el 50% de curvas más profundas (Rango Intercuartílico) y tener en cuenta la función mediana. Luego, para hallar los bigotes primero se encuentran los outliers, es decir, las curvas que se encuentran más allá de \(1.5\) veces el RIQ. Excluyendo estas observaciones, se encuentran los bigotes tomando las líneas verticales de los inliers, como se muestra a continuación en la Figura 10.
Código
p <- 0.5
n <- 267
region <- ceiling(n / 2)
Y <- t(onda_240)
depths <- depths[order(depths$MBD, decreasing = T), ]
BD50 <- Y[depths$idx,][1:region, ]
Min <- apply(BD50, 2, min)
Max <- apply(BD50, 2, max)
RIQ <- Max - Min
lower <- Min - 1.5 * RIQ
upper <- Max + 1.5 * RIQ
isOut <- apply(Y, 1, function(f) any(f < lower | f > upper))
Ygood <- Y[!isOut,]
lowerWhisker <- apply(Ygood, 2, min)
upperWhisker <- apply(Ygood, 2, max)
plot(fd_240, col = "lightgrey", lwd = 1.5,
xlab = "Longitud de Onda de Excitación (nm)",
ylab = "Intensidad de Fluorescencia",)[1] "done"Código
for (k in 1:sum(isOut)) {
lines(x, Y[which(isOut)[k],], col = "red", lty = "dashed")
}
lines(x, eval.fd(x,fd_240_med), col = "black", lwd = 2)
lines(x, Min, col = "springgreen3", lwd = 3)
lines(x, Max, col = "springgreen3", lwd = 3)
polygon(x = c(x, rev(x)),
y = c(Max, rev(Min)),
col = adjustcolor("springgreen", alpha.f = 0.25),
border = NA)
lines(x, lowerWhisker, col = "springgreen4", lwd = 3)
lines(x, upperWhisker, col = "springgreen4", lwd = 3)
Las observaciones atípicas (outliers) están marcadas de color rojo. Estas corresponden a las curvas que se muestran en la Tabla 1.
| Curva | MBD |
|---|---|
| 131 | 0.0997 |
| 10 | 0.0693 |
| 16 | 0.0388 |
| 17 | 0.0343 |
| 14 | 0.0334 |
| 71 | 0.0294 |
Para una mejor visualización gráfica, se muestran únicamente las curvas outliers:
[1] "done"
4.9 Outliers usando el boxplot funcional ajustado
En primer lugar, para hallar el valor de F adecuado, se debe estimar la estructura de covarianza de forma robusta. Para esto, se utiliza el siguiente estimador para la matriz de varianza covarianza:
\[ \operatorname{Cov}(X,Y) = \frac{\sigma_X \sigma_Y}{4} \left[ \operatorname{Var}\!\left(\frac{X}{\sigma_X} + \frac{Y}{\sigma_Y}\right) - \operatorname{Var}\!\left(\frac{X}{\sigma_X} - \frac{Y}{\sigma_Y}\right) \right] \]
donde \(\sigma\) se estima de la siguiente forma:
\[ \hat{\sigma}_Z = Q_n(Z) = d \left\{\, |z_i - z_j| \; ; \; i < j,\; i,j = 1,2,\ldots,n \right\}_{(k)} \] con \(k = \left\lfloor \frac{\binom{n}{2} + 2}{4} \right\rfloor + 1\).
Código
# Estimación de la varianza robusta
var_robusta <- function(z){
n = length(z)
d = 2.2191
k = floor((choose(n,2)+2)/4) + 1
distancia <- c()
for (i in 1:(n-1)){
for (j in (i+1):n){
distancia <- c(distancia, abs(z[i] - z[j]))
}
}
distancia <- sort(distancia, decreasing = FALSE)
var_est <- d*distancia[k]
return(var_est)
}
# Estimación de la matriz de covarianza
p <- 571
var_rob <- numeric(p)
k <- floor((choose(n,2)+2)/4) + 1
X_240 <- wvl[ , , 2]
X_240 <- t(X_240)
for(i in 1:p){
var_rob[i] <- Qn(X_240[, i], k=k)
}
matriz_cov <- matrix(NA, nrow = p, ncol = p)
for (i in 1:p){
matriz_cov[i,i] <- (var_rob[i])^2
for (j in i:p){
if(i != j){
factor1 <- (var_rob[i] * var_rob[j]) / 4
factor2 <- (Qn((X_240[,i] / var_rob[i]) +
(X_240[,j] / var_rob[j]), k = k))^2 -
(Qn((X_240[,i] / var_rob[i]) -
(X_240[,j] / var_rob[j]), k = k))^2
matriz_cov[i,j] <- factor1 * factor2
matriz_cov[j,i] <- matriz_cov[i,j]
}
}
}Sin embargo, la matriz estimada no es definida positiva
Código
sum(eigen(matriz_cov)$values < 0)[1] 275Por lo tanto, se busca una matriz definida positiva cercana a la estimación de la matriz de covarianza
Código
matriz_cov_pd <- as.matrix(nearPD(matriz_cov)$mat)A partir de la matriz definida positiva más cercana a la estimación de la matriz de covarianza, se determina mediante simulación el valor de \(F\) tal que la proporción estimada de observaciones clasificadas como inliers por el boxplot sea aproximadamente del \(99.3\%\).
Código
F_grid <- seq(1.2, 2, by = 0.05)
v_exp <- seq(275, 560, by = 0.5)
n_obs <- 100
n_sim <- 1000
cobertura <- matrix(NA, nrow = n_sim, ncol = length(F_grid))
for (i in 1:n_sim){
sim <- mvrnorm(n_obs, rep(0, 571), matriz_cov_pd)
perfil_sim <- smooth.basis(argvals = v_exp,
y = t(sim),
fdParobj = fd_240)
curvas_eval_sim <- eval.fd(v_exp, perfil_sim$fd)
depth <- MBD(t(curvas_eval_sim))
index <- order(depth, decreasing = TRUE)
m <- ceiling(n_obs * 0.5)
center <- curvas_eval_sim[, index[1:m]]
inf <- apply(center, 1, min)
sup <- apply(center, 1, max)
for (j in seq_along(F_grid)){
f <- F_grid[j]
dist <- f * (sup - inf)
upper <- sup + dist
lower <- inf - dist
outly <- (curvas_eval_sim <= lower) | (curvas_eval_sim >= upper)
outcol <- colSums(outly)
cobertura[i,j] <- (n_obs - sum(outcol > 0)) / n_obs
}
}
cobertura_media <- colMeans(cobertura)
df_cobertura <- data.frame(
F = F_grid,
cobertura = cobertura_media
)Código
df_cobertura <- read.csv("df_cobertura.csv")
df_cobertura F cobertura
1 1.20 0.98149
2 1.25 0.98488
3 1.30 0.98755
4 1.35 0.99000
5 1.40 0.99175
6 1.45 0.99339
7 1.50 0.99474
8 1.55 0.99584
9 1.60 0.99667
10 1.65 0.99737
11 1.70 0.99783
12 1.75 0.99817
13 1.80 0.99852
14 1.85 0.99880
15 1.90 0.99908
16 1.95 0.99923
17 2.00 0.99934Con estos resultados, se toma \(F=1.45\), pues la proporción estimada de observaciones clasificadas como inliers por el boxplot funcional es de \(0.99339\), el cual es el valor más cercano a \(0.993\). El boxplot funcional ajustado se ilustra en la Figura 12.
Código
fda::fbplot(curvas_eval, method ="MBD", factor = 1.45, xlab="Longitud de Onda de Excitación (nm)", ylab = "Intensidad de Fluorescencia", main = "Boxplot Funcional Ajustado")$depth
[1] 0.40437035 0.26123027 0.18343636 0.32432633 0.32846491 0.35339746
[7] 0.30715304 0.25780843 0.22882029 0.06872910 0.24461042 0.33299822
[13] 0.13786754 0.03227213 0.08440709 0.03835688 0.03358344 0.11697882
[19] 0.07967275 0.26690390 0.17040619 0.37435301 0.48580655 0.42541052
[25] 0.48340010 0.48048327 0.43053816 0.48450994 0.45292007 0.38962841
[31] 0.34123301 0.38618344 0.45271905 0.43965647 0.40190664 0.42836499
[37] 0.30815374 0.06176409 0.14966429 0.42629572 0.21928121 0.42943093
[43] 0.43172627 0.48898669 0.17946241 0.23626339 0.32611429 0.16626722
[49] 0.25053020 0.42117203 0.45418274 0.43519324 0.48683659 0.47790229
[55] 0.45822323 0.27219439 0.22025025 0.39730034 0.38990277 0.36115368
[61] 0.27352922 0.30995931 0.15006869 0.12620992 0.32071116 0.47147903
[67] 0.42779793 0.17711031 0.33331933 0.36735713 0.02877967 0.41180664
[73] 0.45923369 0.45352021 0.39912257 0.41281937 0.39764601 0.45269138
[79] 0.49342541 0.34564974 0.24648735 0.46683233 0.41429791 0.38779548
[85] 0.28349105 0.27954437 0.17984995 0.14897093 0.16056124 0.32807865
[91] 0.16674959 0.12864340 0.41591528 0.41648006 0.36831324 0.33109161
[97] 0.49313597 0.20825574 0.46531321 0.44029494 0.48153630 0.42474592
[103] 0.48728509 0.41448990 0.25175194 0.42584777 0.45863330 0.37010455
[109] 0.34903509 0.36035315 0.45676289 0.12544979 0.42237646 0.32248260
[115] 0.44087772 0.43113066 0.33943110 0.30657993 0.26275280 0.36490846
[121] 0.41251252 0.37335941 0.17457386 0.33627961 0.41203601 0.44377803
[127] 0.34553887 0.39084478 0.05388779 0.17578318 0.09842918 0.27720588
[133] 0.24404347 0.37490127 0.34193253 0.19081357 0.07452726 0.47195992
[139] 0.48563384 0.45644257 0.47401572 0.45900964 0.45563712 0.46549844
[145] 0.46245156 0.42638563 0.45620762 0.35535951 0.41172960 0.33563473
[151] 0.35258959 0.43275168 0.37518337 0.37122722 0.43927303 0.45548117
[157] 0.40524677 0.44044294 0.45671815 0.45043511 0.45186512 0.23581638
[163] 0.27308077 0.22105284 0.36931967 0.16753483 0.30838085 0.43158009
[169] 0.46482511 0.36023119 0.44894231 0.14838933 0.35141051 0.23881014
[175] 0.49192113 0.46761747 0.32670422 0.37856827 0.47259350 0.24972420
[181] 0.48435469 0.47819730 0.45757302 0.48329693 0.33375825 0.22054837
[187] 0.15796847 0.12305479 0.26900478 0.33179379 0.27599583 0.36128190
[193] 0.39218479 0.38651959 0.46066069 0.39786024 0.44188750 0.40417214
[199] 0.35269326 0.41199769 0.40232569 0.41318190 0.42842648 0.37811825
[205] 0.22427228 0.38296177 0.36565952 0.25420564 0.31679787 0.37473843
[211] 0.38462634 0.16781535 0.40596212 0.25283481 0.35543541 0.38105767
[217] 0.38168258 0.17084595 0.41755528 0.46087547 0.46702975 0.25459510
[223] 0.32282496 0.15142542 0.14706974 0.27649586 0.30272571 0.43604559
[229] 0.43694465 0.43404616 0.44186037 0.39698076 0.24529258 0.23994376
[235] 0.40976105 0.17032442 0.20507604 0.31814039 0.21654221 0.37782427
[241] 0.26381934 0.29680732 0.39932413 0.45011464 0.42304945 0.37954077
[247] 0.32804645 0.20923780 0.22956188 0.33932975 0.34887471 0.36161499
[253] 0.37665342 0.15975825 0.24850394 0.33988467 0.37805656 0.39227400
[259] 0.37576615 0.34995234 0.32234125 0.28399557 0.32787921 0.28583876
[265] 0.42619585 0.41979578 0.34782054
$outpoint
[1] 10 14 16 17 71 131
$medcurve
[1] 79
4.10 Outliers usando el FOM y el mapa de calor
De acuerdo a lo propuesto por Rousseeuw, Raymaekers, y Hubert (2018) se presenta el concepto de atipicidad.
Stahel–Donoho Outlyingness (SDO)
En un caso univariado, el SDO de un punto \(y\) respecto a una muestra \(Y = \{y_1, y_2, \dots, y_n\}\) se define como:
\[ \begin{equation}\text{SDO}(y; Y) = \frac{|y - \text{med}(Y)|}{\text{MAD}(Y)}\end{equation} \]
Donde el denominador es la Desviación Absoluta de la Mediana (MAD):
\[ \begin{equation}\text{MAD}(Y) = \frac{\text{med}_i \left( |y_i - \text{med}_j(y_j)| \right)}{\Phi^{-1}(0.75)}\end{equation} \]
Una limitación del SDO es que asume implícitamente que los datos no atípicos (inliers) siguen una distribución aproximadamente simétrica. Sin embargo, cuando la distribución de los inliers es asimétrica, el uso del MAD como única medida de escala no logra capturar dicha asimetría. Por ejemplo, si los datos provienen de una distribución sesgada hacia la derecha, algunos inliers ubicados en la cola derecha pueden presentar valores elevados de outlyingness y ser considerados atípicos, mientras que verdaderos outliers ubicados en la cola izquierda podrían no ser detectados adecuadamente.
Directional outlyingness (DO)
Sea \(y_1 \leq y_2 \leq \dots \leq y_n\) una muestra univariada, se contruye dos muestra de tamaño \(h = \left\lfloor \frac{n+1}{2} \right\rfloor\)
Si \(h\) es par entonces se toma \(Y_a = \{y_{h+1}, \dots ,y_n \}\) y \(Y_b = \{y_1, \dots y_h\}\)
Si \(h\) es impar entonces se toma \(Y_a = \{y_{h}, \dots ,y_n \}\) y \(Y_b = \{y_1, \dots y_h\}\)
Inicialmente, se calculan los siguientes estimadores de dispersión.
\[ \begin{equation} s_{o,a}(Y) = \frac{\text{med}(Z_a)}{\Phi^{-1}(0.75)} \quad \text{y} \quad s_{o,b}(Y) = \frac{\text{med}(Z_b)}{\Phi^{-1}(0.75)} \end{equation} \] Donde \(Z_a = Y_a - \text{med}(Y)\) y \(Z_b = Y_b - \text{med}(Y)\).
Luego se calculan las escalas finales \(s_a\) y \(s_b\) utilizando la función \(\rho_c\) de Huber
\[ \begin{equation}s_a(Y) = s_{o,a}(Y) \sqrt{ \frac{1}{2\alpha h} \sum_{z_i \in Z_a} \rho_c \left( \frac{z_i}{s_{o,a}(Y)} \right) }\end{equation} \]
\[ \begin{equation}s_b(Y) = s_{o,b}(Y) \sqrt{ \frac{1}{2\alpha h} \sum_{z_i \in Z_b} \rho_c \left( \frac{z_i}{s_{o,b}(Y)} \right) }\end{equation} \]
\[ \alpha = \int_{0}^{\infty} \rho_c(x) d\Phi(x) \qquad \text{y} \qquad \rho_c(t) = \left( \frac{t}{c} \right)^2 \mathbb{I}_{[-c,c]} + \mathbb{I}_{(-\infty, -c) \cup (c, \infty)} \]
Asi, el DO de un punto \(y\) respecto a una muestra \(Y = \{y_1, y_2, \dots, y_n\}\) se define como:
\[ \begin{equation}\text{DO}(y; Y) = \begin{cases} \frac{y - \text{med}(Y)}{s_a(Y)} & \text{si } y \geq \text{med}(Y) \\\frac{\text{med}(Y) - y}{s_b(Y)} & \text{si } y \leq \text{med}(Y)\end{cases}\end{equation} \]
El MDO se define mediante proyecciones
\[ \text{MDO}(\mathbf{y}; \mathbf{Y}) = \sup_{\mathbf{v} \in \mathbb{R}^d} \text{DO}(\mathbf{y}^T \mathbf{v}; \mathbf{Y}^T \mathbf{v}) \]
Functional Directional Outlyingness
Dado un conjunto de datos funcionales como \(\mathbf{Y} = \{Y_1, Y_2, \dots, Y_n\}\), donde cada \(Y_i\) es una función. Tal que estas funciones se observan en un conjunto discreto de puntos de su dominio \(\{t_1, \dots, t_T\}\).
Para define la medida de alejamiento funcional direccional (fDO) de \(X\) respecto a la muestra \(\mathbf{Y}\) como un promedio ponderado de los alejamientos en cada punto.
\[ \text{fDO}(X; \mathbf{Y}) = \sum_{j=1}^{T} \text{DO}(X(t_j); \mathbf{Y}(t_j)) W(t_j) \]
Donde \(W(\cdot)\) es una función de pesos para la cual:
\[ \sum_{j=1}^{T} W(t_j) = 1 \]
The Functional Outlier Map
El FOM se define como el diagrama de dispersión de que permite identificar de forma simultánea tanto la atipicidad de cada función como la variabilidad de dicho atipicidad a lo largo del dominio: \[(\text{fDO}(Y_i; \mathbf{Y}), \text{vDO}(Y_i; \mathbf{Y}))\]
La ecuación que mide la variabilidad de la atipicidad funcional es:
\[ \begin{equation}\text{vDO}(Y_i; \mathbf{Y}) = \frac{\text{stdev}_j (\text{DO}(Y_i(t_j); \mathbf{Y}(t_j)))}{1 + \text{fDO}(Y_i; \mathbf{Y})}\end{equation} \]
Para señalar valores atípicos se define combined functional outlyingness (CFO de una función Yi como
\[ \begin{equation}\begin{array}{ll}\text{CFO}_i & = \text{CFO}(Y_i; \mathbf{Y}) \\\\& = \sqrt{(\text{fDO}_i / \text{med}(\text{fDO}))^2 + (\text{vDO}_i / \text{med}(\text{vDO}))^2}\end{array}\end{equation} \]
Para determinar un cutoff para \(\text{CFO}\), para cada \(i = 1, \dots, n\), se calcula: \[ LCFO_i = \log(0.1 + \text{CFO}_i) \]Se considera la función \(Y_i\) como atípica si se cumple que:
\[ \begin{equation}\frac{LCFO_i - \text{med}(LCFO)}{\text{MAD}(LCFO)} > \Phi^{-1}(0.995)\end{equation} \]
Usando la funcion fOutl del paquete mrfDepth se calculan las atipicidades direccionales univariadas y con la funcion fom el Functional Outlier Map.
Código
matriz_240 <- t(onda_240)
array_240 <- array(matriz_240, dim = c(268, 571,1))
DO_240 <- fOutl(x = aperm(array_240, c(2, 1, 3)), alpha = 0,type = "fDO", distOptions = list(rmZeroes = TRUE,
maxRatio = 3),
diagnostic = TRUE)
fom(DO_240, cutoff = TRUE)
Como se puede observar en Figura 13, seis puntos son identificados como atípicos, cuando el valor del fDO es alto y la medida de variación de la atipicidad (vFO) es pequeña, esto sugiere que la función presenta un nivel elevado de atipicidad que se mantiene de manera relativamente constante a lo largo del dominio.
Por otro lado, los puntos que se encuentran dentro de la región representada por una elipse, corresponden a funciones que, si bien pueden presentar valores ligeramente mayores o menores de atipicidad, no muestran una variabilidad elevada en esta medida. Esto sugiere que dichas funciones no se apartan significativamente del comportamiento general del conjunto, o bien que las posibles desviaciones corresponden a atipicidades de carácter más local.
Finalmente, dado que la implementación del FOM proporcionada por el paquete utilizado no permite identificar directamente los índices de las observaciones atípicas, se procederá a construir el FOM desde cero para identificar las observaciones que pueden considerarse atípicas dentro del conjunto de datos.
Código
fDO = function(DO,weights=NULL){
if(is.null(weights)){
weights=apply(DO,MARGIN = 2:(length(dim(DO))),FUN = function(y)
(prod(is.na(y))==0))
}
weights = weights/sum(weights,na.rm = TRUE)
fDOvalues = apply(DO,1,FUN= function(y) sum(weights*y,na.rm=TRUE))
return(fDOvalues)
}
fdo240 <- fDO(DO_240$crossDistsX)
desv240 <- apply(DO_240$crossDistsX, 1, sd)
vDO240 <- desv240/(1+fdo240)
CFO240 <- sqrt(( fdo240/median(fdo240))^2 + (vDO240/median(vDO240))^2)
LCFO240 <- log(0.1 + CFO240)
z_LCFO240 <- (LCFO240 - median(LCFO240)) / mad(LCFO240)
umbral <- qnorm(0.995)
es_outlier_cfo <- z_LCFO240 > umbral
indices_a_marcar <- which(es_outlier_cfo)
x_grid <- seq(0, max(fdo240) * 1.2, length.out = 200)
y_grid <- seq(0, max(vDO240) * 1.2, length.out = 200)
grid <- expand.grid(x = x_grid, y = y_grid)
norm_grid_x <- grid$x / median(fdo240)
norm_grid_y <- grid$y / median(vDO240)
CFO_grid <- sqrt(norm_grid_x^2 + norm_grid_y^2)
LCFO_grid <- log(0.1 + CFO_grid)
z_LCFO_grid <- (LCFO_grid - median(LCFO240)) / mad(LCFO240)
z_matrix <- matrix(z_LCFO_grid, nrow = length(x_grid))
contorno <- contourLines(x_grid, y_grid, z_matrix, levels = umbral)
colores_puntos <- ifelse(es_outlier_cfo, "red", "black")
plot(fdo240, vDO240,
col = colores_puntos,
pch = 19,
cex = 0.8,
xlab = "fDO",
ylab = "sd(fDO) / (1+fDO)",
xlim = c(0, max(fdo240) * 1.2),
ylim = c(0, max(vDO240) * 1.2),
axes = TRUE
)
if(length(contorno) > 0) {
lines(contorno[[1]]$x, contorno[[1]]$y, lty = 2, lwd = 2, col = "black")
}
text(fdo240[indices_a_marcar], vDO240[indices_a_marcar],
labels = indices_a_marcar,
cex = 0.5, col = "red", pos = 3, offset = 0.5)
Como se puede observar en Figura 14, se conserva la misma estructura obtenida previamente al utilizar el paquete, lo que confirma la consistencia de los resultados. A partir de esta construcción, es posible identificar explícitamente las observaciones atípicas, las cuales corresponden a las funciones 71, 10, 13, 14, 16, 129 y 131.
Finalmente, se construira el mapa de calor de los DO. Este análisis permitirá examinar cómo se distribuye la atipicidad a lo largo del dominio. En particular, dado que se cuenta con 571 mediciones de longitud de onda, el mapa de calor permitirá identificar las regiones del dominio donde los valores de DO son mayores, así como las observaciones que concentran estos niveles elevados de atipicidad direccional.
Código
DO <- DO_240$crossDistsX
rownames(DO) <- c(1:268)
colnames(DO) <- x
cap <- quantile(DO, 0.99)
DO[DO > cap] <- cap
df <- as.data.frame(as.table(DO))
colnames(df) <- c("Observacion", "Longitud", "DO")
df$Longitud <- as.numeric(as.character(df$Longitud))
df$Observacion <- as.numeric(as.character(df$Observacion))
ggplot(df, aes(x = Longitud, y = Observacion, fill = DO)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "darkred") +
scale_x_continuous(
breaks = seq(min(df$Longitud), max(df$Longitud), by = 50)
) +
scale_y_continuous(
breaks = seq(min(df$Observacion), max(df$Observacion), by = 50)
) +
theme_classic()
En Figura 15 se observan varias lineas horizontales con valores elevados, lo que indica la presencia de observaciones que mantienen niveles altos de atipicidad direccional a lo largo de gran parte del dominio de longitudes de onda. Este patrón sugiere que la atipicidad identificada no se concentra únicamente en regiones específicas, sino que en algunos casos corresponde a atipicidades de carácter más global, donde ciertas funciones se mantienen alejadas del comportamiento central en múltiples longitudes de onda. Asimismo, se identifican zonas con valores de DO relativamente altos en distintos intervalos de observaciones, lo cual podría corresponder a algunas de las funciones previamente detectadas como atípicas mediante el FOM, tales como las observaciones 14, 16, 71, 129 y 131.
5 Profundidad Multivariada
Para el cálculo de la profundidad multivariada para los siete procesos funcionales, se debe tener en cuenta que se tiene una muestra de tamaño \(n=268\), \(\{\mathbf{Y}_1, \mathbf{Y}_2, \ldots, \mathbf{Y}_{268} \}\) donde \(U=\mathcal{T}=[275, 560]\) y \(\boldsymbol{{Y}}_i \in \mathcal{C}(U)^7\) con \(i=1,\ldots, 268\). El cálculo de la profundidad se realiza mediante la MFHD (multivariate functional halfspace depth).
MFHD corresponde a la Multivariate Functional Halfspace Depth. La cual está dada por:
\[MFHD_N(\boldsymbol{\chi})=\sum_{j=1}^T HD(\boldsymbol{\chi}(t_j); F_{\boldsymbol{\mathcal{Y}}(t_j), N})W_j\]
donde la profundidad \(HD\) muestral está dada por:
\[HD(X(t); F_{\mathcal{Y}(t),N}) =\frac{1}{N}\min_{u \in\mathbb{R}^K,\ \|u\|=1}\#\left\{ Y_n(t),\ n = 1,\ldots,N :u'Y_n(t) \ge u'X(t) \right\}\]
Esta medida de profundidad consiste en dividir el espacio en dos partes mediante un hiperplano que pasa por la observación \(\boldsymbol{X}(t)\). Para cada dirección \(u\), los datos se proyectan sobre dicha dirección y se calcula la proporción de observaciones \(Y_n(t)\) que satisfacen \(u'Y_n(t) \ge u'X(t)\). La profundidad se define como la mínima proporción obtenida al considerar todas las direcciones posibles. La idea intuitiva de esta medida de profundidad es medir cuántos datos se encuentran alrededor del valor \(\boldsymbol{X}(t)\), lo que permite determinar qué tan central es esta observación respecto al conjunto de datos.
En conclusión, la MFHD consiste en calcular la HD de un vector de funciones \(\boldsymbol{\chi}\) en la muestra punto a punto a lo largo de su dominio, con respecto a los demás vectores de curvas, y luego obtener un promedio ponderado de estas profundidades.
Código
fitted.curves <- function(fd, t){
fd_obj <- fd
grid <- t
Yt <- eval.fd(grid, fd_obj)
Y <- t(Yt)
return(list(Y = Y, Yt = Yt))
}Código
YY <- array(NA, dim = c(1000, 267, 7)) #t by n by p array
# t time points
# n functional observations
# p number of measurements for every functional observation
grid <- seq(275, 560, length = 1000)
resultados_suavizados <- resultados_suavizados[order(as.numeric(names(resultados_suavizados)))]
names(resultados_suavizados) <- 1:7
for (k in 1:7) {
YY[, , k] <- fitted.curves(resultados_suavizados[[k]]$fd, t = grid)$Yt
}Código
MFHD <- mfd(YY, type = "hdepth")
Mdepths <- data.frame(idx = 1:length(MFHD$MFDdepthX),
depth = MFHD$MFDdepthX)
Mdepths <- Mdepths[order(Mdepths$depth, decreasing = TRUE),]| Curva | MBD |
|---|---|
| 34 | 0.0568 |
| 184 | 0.0566 |
| 195 | 0.0548 |
| 33 | 0.0538 |
| 75 | 0.0524 |
| 111 | 0.0523 |
| 192 | 0.0499 |
| 204 | 0.0489 |
| 28 | 0.0483 |
| 27 | 0.0479 |
Orden de los datos funcionales multivariados con la MFHD
5.1 Función Mediana Multivariada
Para el caso multivariado, la función mediana corresponde al vector \(\mathbf{Y}_i\) que maximiza la MFHD. En este caso, dicho vector es \(\mathbf{Y}_{34}\), y se muestra en la Figura 16.
5.2 Observaciones Outliers
Sea \(\Psi=\{\psi_1,\ldots,\psi_n\}\) un conjunto de datos funcionales multivariados, donde \(\psi_i = (\psi_{i1},\ldots,\psi_{ip}), \quad i=1,\ldots,n\) y \(\psi_{ij}\) está definida en \([a,b]\) para \(j=1,\ldots,p\).
Sea \(\chi = (\chi_1, \ldots ,\chi_p)'\) donde \(\chi_j \subset [a,b]\), \(j=1,\ldots,p\).
Sea \(\{t_1,\ldots,t_N\}\) un conjunto discreto de puntos tales que \(t_j \subset [a,b]\).
La atipicidad direccional (DO) de la función multivariada \(\chi\) con respecto a \(\Psi\) está dada por
\[ \begin{equation} \text{mFDO}(\chi,\Psi) = \sum_{j=1}^{N} \text{MDO}\big(\chi(t_j);\Psi(t_j)\big)\,\omega(t_j) \end{equation} \]
De esta forma, igual que en el caso univariado, se puede graficar el FOM y establecer un punto de corte para la detección de valores atípicos.
Código
# Directional Outlyiness
f_out <- fOutl(x = wvl, type = "fDO", diagnostic = TRUE)
fom_X <- fom(f_out, cutoff = TRUE)Código
indices_outliers_fom <- function(fOutlResult){
AOValues <- fOutlResult$crossDistsX
weights <- fOutlResult$weights
fAO <- AOValues %*% weights
sdAO <- apply(AOValues, 1, function(y){
sqrt(sum(weights * (y - sum(weights * y, na.rm=TRUE))^2, na.rm=TRUE) /
(1 - 1/length(weights)))
})
DispMeasure <- sdAO / (1 + fAO)
CAO <- log(0.1 + sqrt((fAO/median(fAO))^2 +
(DispMeasure/median(DispMeasure))^2))
Fence <- qnorm(0.995) * mad(CAO) + median(CAO)
theta <- seq(0, pi/2, length = (100))
FenceData <- matrix(0, nrow = length(theta), ncol = 2)
colnames(FenceData) <- c("x", "y")
FenceData <- data.frame(FenceData)
FenceData$x <- median(fAO) * (exp(Fence) - 0.1) * cos(theta)
FenceData$y <- median(DispMeasure) * (exp(Fence) - 0.1) * sin(theta)
PlotData <- ifelse(CAO > Fence, "red", "black")
outliers <- which(CAO > Fence)
plot(fAO, DispMeasure,
col = PlotData,
pch = 16,
xlab = "MfDO", ylab = "MvDO")
lines(FenceData$x, FenceData$y, lty = 2, lwd = 2)
text(fAO[outliers], DispMeasure[outliers],
labels = outliers,
col = "red",
pos = 3, cex = 0.5)
}
indices_outliers_fom(f_out)
6 Replique la sección 4. Application to Image Data del artículo de Rousseeuw, P. J.,Raymaekers, J., & Hubert, M. (2018).
Se replicara la sección 4 del artículo Rousseeuw, Raymaekers, y Hubert (2018), donde se trabaja una aplicación de la metodología de Directional Outlyingness con imágenes.
Las imágenes son funciones sobre un dominio bivariado. En la práctica, el dominio es una cuadrícula de puntos discretos, por ejemplo, los píxeles horizontales y verticales de una imagen. Por lo que es conveniente usar dos índices \(j = 1, \dots, J\) y \(k = 1, \dots, K\), uno para cada dimensión de la cuadrícula. Los valores de la función pueden ser univariados, como intensidades de gris, pero también pueden ser multivariados, por ejemplo, las intensidades de rojo, verde y azul (RGB).
En general, se escribira un conjunto de datos de imágenes como una muestra \(\boldsymbol{Y} = \{Y_1, Y_2, \dots, Y_n\}\) donde cada \(Y_i\) es una función de \(\{(j, k); j = 1, \dots, J \text{ y } k = 1, \dots, K\}\) a \(\mathbb{R}^d\).
Las nociones de \(\text{fDO}\) y \(\text{vDO}\) desarrolladas anteriormente seran extendidas, así;
\[\begin{equation} \text{fDO}(Y_i; \boldsymbol{Y}) = \sum_{j=1}^{J} \sum_{k=1}^{K} \text{DO}(Y_i(j, k); \boldsymbol{Y}(j, k))W_{jk}, \end{equation}\]
donde los pesos \(W_{jk}\) deben satisfacer \(\sum_{j=1}^{J} \sum_{k=1}^{K} W_{jk} = 1\), y
\[\begin{equation} \text{vDO}(Y_i; \boldsymbol{Y}) = \frac{\text{stdev}_{j,k}(\text{DO}(Y_i(j, k); \boldsymbol{Y}(j, k)))}{1 + \text{fDO}(Y_i; \boldsymbol{Y})}, \end{equation}\]
Se estudiara un dataset que contiene imágenes de resonancia magnética cerebral de 416 sujetos entre 18 y 96 años. Las imágenes tienen \(176\) por \(208\) píxeles, con valores de escala de grises entre \(0\) y \(255\). En conjunto, tenemos por tanto 416 imágenes observadas \(Y_i\) que contienen valores de intensidad univariados \(Y_i(j, k)\), donde \(j = 1, \dots, J = 176\) y \(k = 1, \dots, K = 208\).
Código
load("DO_MRI_data.RData")Usando la función fOutl del paquete mrfDepth se calculara el DO punto por punto.
Código
dims <- dim(MRI)
dim(MRI) <- c(dims[1], prod(dims[2:3]), dims[4])
DO_MRI_compwise <- mrfDepth::fOutl(x = aperm(MRI, c(2, 1, 3)), type = "fDO", distOptions = list(maxRatio = 2, type = "compWise"), diagnostic = TRUE)También se grafica el FOM en este caso con el fin de identificar los vectores de imágenes que presentan un comportamiento atípico.
Código
indices_outliers_fom <- function(fOutlResult){
AOValues <- fOutlResult$crossDistsX
weights <- fOutlResult$weights
fAO <- AOValues %*% weights
sdAO <- apply(AOValues, 1, function(y){
sqrt(sum(weights * (y - sum(weights * y, na.rm=TRUE))^2, na.rm=TRUE) /
(1 - 1/length(weights)))
})
DispMeasure <- sdAO / (1 + fAO)
CAO <- log(0.1 + sqrt((fAO/median(fAO))^2 +
(DispMeasure/median(DispMeasure))^2))
Fence <- qnorm(0.995) * mad(CAO) + median(CAO)
theta <- seq(0, pi/2, length = (100))
FenceData <- matrix(0, nrow = length(theta), ncol = 2)
colnames(FenceData) <- c("x", "y")
FenceData <- data.frame(FenceData)
FenceData$x <- median(fAO) * (exp(Fence) - 0.1) * cos(theta)
FenceData$y <- median(DispMeasure) * (exp(Fence) - 0.1) * sin(theta)
PlotData <- ifelse(CAO > Fence, "red", "black")
outliers <- which(CAO > Fence)
plot(fAO, DispMeasure,
col = PlotData,
pch = 16,
xlab = "MfDO", ylab = "MvDO")
lines(FenceData$x, FenceData$y, lty = 2, lwd = 2)
text(fAO[outliers], DispMeasure[outliers],
labels = outliers,
col = "red",
pos = 3, cex = 0.5)
}
indices_outliers_fom(DO_MRI_compwise)
7 Referencias
López-Pintado, Sara, y Juan Romo. 2009. «On the Concept of Depth for Functional Data». Journal of the American Statistical Association 104 (486): 718-34. https://doi.org/10.1198/jasa.2009.0108.
Ramsay, J. O., Giles Hooker, y Spencer Graves. 2009. Functional Data Analysis with R and MATLAB. Springer.
Rousseeuw, Peter J., Jakob Raymaekers, y Mia Hubert. 2018. «A Measure of Directional Outlyingness With Applications to Image Data and Video». Journal of Computational and Graphical Statistics 27 (2): 345-59. https://doi.org/10.1080/10618600.2017.1366912.