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)
library(matrixcalc)
library(ccaPP)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)\).
sugar <- readMat("C:\\Users\\sofia\\OneDrive\\Documentos\\2026-1\\FDA\\data.mat")
sugar_fluo <- sugar$X
dim(sugar_fluo) [1] 268 3997t <-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)
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 la función smooth.pos del paquete fda para garantizar que el ajuste no de valores negativos y \(K=40\).
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 la suma de cuadrados penalizada
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]))ondas_lista <- list(
"340" = onda_340, "325" = onda_325, "305" = onda_305,
"290" = onda_290, "255" = onda_255, "240" = onda_240,
"230" = onda_230
)
lambdas <- 10^seq(-3, 4, 0.25)
BSpl <- create.bspline.basis(rangeval = c(275, 560), nbasis = 40)
ncurve <- dim(wvl)[2]
for (nombre in names(ondas_lista)) {
data_actual <- ondas_lista[[nombre]]
gcv_values <- numeric(length(lambdas))
for(i in 1:length(lambdas)){
Wfd0 <- fd(matrix(0, BSpl$nbasis, ncurve), BSpl)
WfdPar <- fdPar(Wfd0, lambda = lambdas[i])
Fsugar_tmp <- smooth.pos(x, data_actual, WfdPar)
fs <- NULL
for (n in seq_len(ncurve)) {
fs <- c(fs, Fsugar_tmp$Flist[[n]]$f)
}
fsave[i] <- sum(fs)
}
mejor_lambda[nombre] <- lambdas[which.min(fsave)]
cat("Mejor lambda:", mejor_lambda, "\n")
}| Parámetro | Valor Seleccionado |
|---|---|
| \(\lambda_{1}\) | 0.0056 |
| \(\lambda_{2}\) | 0.0010 |
| \(\lambda_{3}\) | 0.0010 |
| \(\lambda_{4}\) | 0.0010 |
| \(\lambda_{5}\) | 1.0000 |
| \(\lambda_{6}\) | 0.0018 |
| \(\lambda_{7}\) | 0.0018 |
Wfd0 <- fd(matrix(0, BSpl$nbasis, ncurve), BSpl)
grid <- seq(min(x), max(x), length = 1000)
fit <- list()
fits <- array(NA, dim = c(length(grid), ncurve, 7))
for (i in 1:7) {
WfdPar <- fdPar(Wfd0, lambda = Lambdas[i])
fit[[i]] <- smooth.pos(x, wvl[, , i], WfdPar, dbglev = 1)
Wfd <- fit[[i]]$Wfdobj
fits[ , , i] <- exp(eval.fd(grid, Wfd))
}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).
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} \]
grid <- seq(min(t), max(t), length = 1000)
fd.240 <- t(fits[ , , 2])
Ybar <- colMeans(fd.240)
matplot(grid, fits[ , , 2],
type = "l",
col = "#CCCCFF",
lty = 1,
xlab = "Longitud de Onda de Excitación (nm)",
ylab = "Intensidad de Fluorescencia")
lines(grid, Ybar, col = "#000033", lwd = 1.5)
legend("topright",
legend = bquote("Función Media " ~ bar(x)(t)),
col = "#000033",
lwd = 1.5,
bty = "n",
cex = 0.9)
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.
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.
fd.240 <- t(fits[ , , 2])
grid <- seq(min(t), max(t), length = 1000)
n <- nrow(fd.240)
alpha <- 0.1
upto <- n - floor(n * alpha)
mbd <- MBD(fd.240, plotting = FALSE)
depths <- data.frame(idx = seq_along(mbd$MBD),
MBD = as.numeric(mbd$MBD))
depths <- depths[order(depths$MBD, decreasing = TRUE), ]
ind <- depths$idx[1:upto]
Ytrimmed <- fd.240[depths$idx,][1:upto,]
Ytrimmed_mean <- colMeans(Ytrimmed)| Índice Curva | MBD |
|---|---|
| 79 | 0.4964 |
| 97 | 0.4922 |
| 44 | 0.4918 |
| 175 | 0.4898 |
| 103 | 0.4880 |
| 53 | 0.4876 |
| 139 | 0.4862 |
| 23 | 0.4858 |
| 181 | 0.4855 |
| 25 | 0.4853 |
Orden de las curvas basadas en la MBD
Con base en estos resultados, se presenta la media recortada al 10% en la Figura 3.
matplot(grid, fits[ , , 2],
type = "l",
col = "#CCCCFF",
lwd = 1.5,
xlab = "Longitud de Onda de Excitación (nm)",
ylab = "Intensidad de Fluorescencia")
lines(grid, Ytrimmed_mean, col = "seagreen3", lwd = 2)
lines(grid, Ybar, col = "#000033", lwd = 2, lty = "solid")
legend("topright",
legend = c(bquote("Función Media " ~ bar(x)(t)),
bquote("Función Media Recortada " ~ bar(x)[0.1](t))),
col = c("#000033",
"seagreen3"),
lty = "solid",
lwd = 2,
bty = "n",
cex = 0.9)
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.
Ybar <- colMeans(fd.240)
centered <- sweep(fd.240, 2, Ybar, "-")
centered2 <- centered^2
var_f <- colSums(centered2) / (nrow(fd.240) - 1)
plot(grid, var_f, type="l", lwd=1,
xlab="Longitud de Onda de Excitación (nm)", ylab="Varianza", col="red", main = "Función Varianza") 
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)]\)
scov_mat <- t(centered) %*% centered / (nrow(fd.240) - 1)
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 = ~grid, y = ~grid, z = ~scov_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"
)
)
)
p