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.
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 que los espectros de emisión no presentan un comportamiento periódico. En este contexto, los B-splines son una base adecuada.
#| code-fold: false
#| echo: true
file <- read.mat("data.mat")
data <- file$X
t <- file$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(data[,start:end])
}
nms <- c("340nm", "325nm", "305nm", "290nm", "255nm", "240nm", "230nm")
nms <- rev(nms)
dimnames(wvl) <- list(NULL, NULL, nms)Basis <- create.bspline.basis(rangeval = c(275, 560),
norder = 6, nbasis = 40)
fit <- list()
grid <- seq(min(t), max(t), length = 1000)
fits <- array(dim = c(length(grid), dim(wvl)[2:3]))
for (k in 1:7) {
fdParobj <- fdPar(fdobj = Basis, Lfdobj = 2, lambda = Lambdas[k])
fit[[k]] <- smooth.basis(t, wvl[, , k], fdParobj)
fdobj <- fit[[k]]$fd
fits[ , , k] <- eval.fd(evalarg = grid, fdobj = fdobj)
}
Para poder realizar FPCA, el operador de covarianza debe ser un operador de Hilbert–Schmidt, simétrico y definido positivo. Para que lo primero se cumpla, el proceso aleatorio del cual se obtienen las realizaciones debe ser cuadrado–integrable, \(\mathrm{E}[\lVert \chi \rVert^2]<\infty\). Las otras dos propiedades se tienen por la simetría y positividad del operador de covarianza \(c(s,t) = \operatorname{Cov}(\chi(s), \chi(t))\).
Otro supuesto que se hace es que la media sea constante. Para esto, se juzga el sistema de hipótesis \[ H_0: \mathrm{E}[\chi_1] = \cdots = \mathrm{E}[\chi_N] \]
Esto se realiza mediante la función fchange() del paquete fChange, donde es posible aplicar la metodología propuesta por (Berkes et al. 2009).
X.dfts <- dfts(fits[ , , 2])
N <- ncol(X.dfts$data)
ChangeDetectionSeg <- fchange(X.dfts,
method = "mean",
statistic = "Mn",
type = "segmentation",
critical = "simulation")A continuación, en la tabla Tabla 1 se muestran la estimaciones de los dos puntos de cambio \(\hat{\tau}_1,\hat{\tau}_2\) junto con sus p–valores.
ChangeDetectionSeg <- read.csv("ChangeDetectionSeg.csv", check.names = FALSE)
colnames(ChangeDetectionSeg) <- c("Punto de Cambio", "p-valor")
ChangeDetectionSeg %>%
kable(format = "html",
digits = 4,
align = "c",
escape = FALSE)| Punto de Cambio | p-valor |
|---|---|
| 134 | 0.000 |
| 184 | 0.004 |
Los resultados proporcionan evidencia para rechazar la hipótesis nula de media funcional constante, \(H_0: \mathrm{E}[\chi_1] = \cdots = \mathrm{E}[\chi_N]\), a favor de la existencia de cambios en la media de las curvas espectrométricas. Se identifican dos puntos de cambio significativos, lo que sugiere la presencia de tres regímenes con distinta estructura de media funcional. En la Figura 1 se muestran las funciones media para cada segmento de curvas y en la se muestra el número de observaciones en cada uno.
ns <- c(cp$Grupo_1$n, cp$Grupo_2$n, cp$Grupo_3$n)
rangoss <- c(cp$Grupo_1$rango, cp$Grupo_2$rango, cp$Grupo_3$rango)
cptable <- data.frame(Segmento = c(1, 2, 3),
n = ns,
Rango = rangoss)
cptable %>%
kable(format = "pandoc",
digits = 4,
align = "c",
escape = FALSE)| Segmento | n | Rango |
|---|---|---|
| 1 | 134 | 1 a 134 |
| 2 | 50 | 135 a 184 |
| 3 | 84 | 185 a 268 |
Los resultados de Tabla 1 proporcionan evidencia para rechazar la hipótesis nula de media funcional constante, a favor de la existencia de cambios en la media de las curvas espectrométricas. Se identifican dos puntos de cambio significativos, lo que sugiere la presencia de tres conjuntos de observaciones con distinta estructura de media funcional.
Por lo tanto, a cada observación se le resta la media estimada correspondiente al segmento al que pertenece, con el fin de obtener las curvas centradas. Este proceso se ilustra en la Figura 2.
A continuación se realiza el análisis de componentes principales para cada uno de los grupos.
X.G1 <- cp$Grupo_1$fd_original
pcaObj.G1 <- pca.fd(X.G1, nharm = 3, centerfns = TRUE)
eigenfn.G1 <- eval.fd(grid, pcaObj.G1$harmonics)
eigenval.G1 <- pcaObj.G1$values[1:3]X.G2 <- cp$Grupo_2$fd_original
pcaObj.G2 <- pca.fd(X.G2, nharm = 3, centerfns = TRUE)
eigenfn.G2 <- eval.fd(grid, pcaObj.G2$harmonics)
eigenval.G2 <- pcaObj.G2$values[1:3]X.G3 <- cp$Grupo_3$fd_original
pcaObj.G3 <- pca.fd(X.G3, nharm = 3, centerfns = TRUE)
eigenfn.G3 <- eval.fd(grid, pcaObj.G3$harmonics)
eigenval.G3 <- pcaObj.G3$values[1:3]X <- t(fits[, , 2])
X.G1 <- t(eval.fd(grid, cp$Grupo_1$fd_original))
X.G2 <- t(eval.fd(grid, cp$Grupo_2$fd_original))
X.G3 <- t(eval.fd(grid, cp$Grupo_3$fd_original))
X.fd <- fit[[2]]
grid <- seq(min(t), max(t), length = 1000)
mean <- colMeans(X)
Xc <- sweep(X, 2, mean, FUN = "-")
mean.Xc <- colMeans(Xc)pcaObj <- pca.fd(X.fd$fd, nharm = 3, centerfns = TRUE)
eigenfn <- eval.fd(grid, pcaObj$harmonics)
eigenval <- pcaObj$values[1:3]Luego de haber centrado los procesos, se encuentran las estimaciones de las funciones propias asociadas al operador de covarianza del proceso funcional para cada grupo.
Las tres primeras funciones propias se pueden visualizar en la Figura 3.
[1] "done"
Con la función plot.pca.fd es posible interpretarlas de una forma más sencilla, pues se muestra cómo cada componente principal afecta a la media del proceso. Esto se muestra en la Figura 4, Figura 5 y Figura 6.
[1] "done"
[1] "done"
En el caso de la primera función propia, esta está mostrando un aumento constante en la media. Debido a su alto porcentaje de variabilidad explicada, se podría decir que la mayor fuente de variación es la cantidad de luz que se emite (magnitud). Además, se puede ver que hay más varibilidad cuando las muestras se someten a longitudes de onda entre los 350nm y 500nm aproximadamente.
La segunda función propia muestra un contraste entre la exposición a las longitudes de onda menores a 380nm aproximadamente y la exposición a las longitudes de onda mayores a 380nm, mostrando una variación mayor en los intervalos [310, 370] y [400, 470] nm.
Por último, la tercera función propia tienen un patrón cuadrático y presenta mayor variabilidad hacia las longitudes de onda medias.
eL porcentaje de variabilidad se muestra a través de un scree plot en la Figura 15.
De esta forma, es posible ver que para los tres grupos, las dos primeras componentes principales ya explican más del 90% de variabilidad total.
La descomposición de Karhunen-Loève consiste en expresar cada observación funcional como \[ X(t) = \mu(t) + \sum_{k=1}^\infty \xi_k \nu_k (t) \] donde \(\nu_k\) corresponde a la \(k\)–ésima función propia y \(\xi_k\) es el \(k\)–ésimo score, \(\xi_k=\langle X, \nu_k \rangle\).
En este caso, la base óptima corresponde al conjunto de las dos primeras funciones propias pues ellas explican más del 90% de la variabilidad, aproximadamente. Por lo tanto, para cada observación funcional, se tiene su aproximación \[ X_{i,g}(t) \approx \hat \mu_g(t) + \sum_{k=1}^{K_g}\hat {\xi}_{ikg} \hat \nu _ {kg}(t) \] con \(i = 1,\ldots,268\), \(g=1,2,3\) y \(K_g=2\).
A continuación, se expresan las tres primeras observaciones muestrales con la descomposición de Karhunen-Loève. En la @ se pueden ver las observaciones junto con su aproximación.
Scores <- matrix(pcaObj$scores[,1:2], nrow = 268, ncol = 2)
eigenfns <- eval.fd(evalarg = grid, fdobj = pcaObj$harmonics)[,1:2]
mu <- eval.fd(evalarg = grid, fdobj = pcaObj$meanfd)
X.approx <- Scores %*% t(eigenfns)
X.approx <- sweep(X.approx, 2, mu, "+") K1 <- 3
Scores.G1 <- matrix(pcaObj.G1$scores[,1:K1], nrow = nrow(X.G1), ncol = K1)
eigenfns.G1 <- eval.fd(evalarg = grid, fdobj = pcaObj.G1$harmonics)[,1:K1]
mu.G1 <- eval.fd(evalarg = grid, fdobj = pcaObj.G1$meanfd)
X.approx.G1 <- Scores.G1 %*% t(eigenfns.G1)
X.approx.G1 <- sweep(X.approx.G1, 2, mu.G1, "+") 
K2 <- 2
Scores.G2 <- matrix(pcaObj.G2$scores[,1:K2], nrow = nrow(X.G2), ncol = K2)
eigenfns.G2 <- eval.fd(evalarg = grid, fdobj = pcaObj.G2$harmonics)[,1:K2]
mu.G2 <- eval.fd(evalarg = grid, fdobj = pcaObj.G2$meanfd)
X.approx.G2 <- Scores.G2 %*% t(eigenfns.G2)
X.approx.G2 <- sweep(X.approx.G2, 2, mu.G2, "+") 
K3 <- 2
Scores.G3 <- matrix(pcaObj.G3$scores[,1:K3], nrow = nrow(X.G3), ncol = K3)
eigenfns.G3 <- eval.fd(evalarg = grid, fdobj = pcaObj.G3$harmonics)[,1:K3]
mu.G3 <- eval.fd(evalarg = grid, fdobj = pcaObj.G3$meanfd)
X.approx.G3 <- Scores.G3 %*% t(eigenfns.G3)
X.approx.G3 <- sweep(X.approx.G3, 2, mu.G3, "+") 
A continuación, en la Tabla 3 se muestra el error cuadrático medio de predicción para cada observación.
| MPSE | |
|---|---|
| Primer Grupo | 81.3371 |
| Segundo Grupo | 770.4291 |
| Tercer Grupo | 294.1825 |
A continuación se presentan las realizaciones del primer score en la ?@fig-primer-score-uni.
Es posible ver que los scores tienen una distribución asimétrica y también presenta algunos outliers por valores muy altos.
A continuación, se muestra el gráfico de dispersión de las dos primeras componentes principales, que explican el 96.33 de la variabilidad, aproximadamente.
Es posible ver que la mayoría de puntos se concentran alrededor de los valores negativos cercanos al cero para el primer score y alrededor del cero para el segundo score. Asimismo, debido al gran porcentaje de variabilidad que explica la primera componente se puede ver que la dispersión es mayor a lo largo de su eje.
En la ?@fig-covarianza-funciones se presenta la estimación de la función de covarianza y las funciones propias.
Como es posible visualizar, en los tres casos la forma de la primera función propia es similar a las que se pueden ver en la función de covarianza debido a su alto porcentaje de variabilidad explicada.
Por el teorema de Mercer se tiene que el kernel \(\psi(s,t)\) usado en el operador integral \(\Psi (x)(t)\) se puede representar en términos de los valores propios y funciones propias del operador. En este caso, se tiene el operador de covarianza \(\mathcal{C}_\chi(x)(t))\) con la función de covarianza como Kernel. Luego,
\[c(s,t) = \sum_{j=1}^\infty \lambda_j \phi(s) \phi(t)\]
En la ?@fig-covarianza-mercer se puede ver la función de covarianza junto con su aproximación y en la ?@fig-covarianza-mercer-diff se ven sus residuales.
Una carta de control es una herramienta que permite monitorear si un proceso se mantiene bajo control o si presenta cambios anómalos a lo largo del tiempo.
Para construir una carta de control a partir del análisis de componentes principales:
Para el caso funcional, sea \(\chi\) una función aleatoria definida en \(L^2\), la cual puede representarse mediante la descomposición de Karhunen–Loève como:
\[\begin{equation} \chi(t) = \mu(t) + \sum_{k=1}^{\infty} \xi_k \, \phi_k(t) \end{equation}\]
Donde:
Para realizar FPCA, se utiliza una aproximación truncada de esta expansión. Dado que el operador de covarianza posee infinitos valores y funciones propias, se selecciona un número finito \(K\) de componentes principales, de manera que se capture una proporción suficientemente alta de la varianza total.
Así, se obtiene la aproximación:
\[\begin{equation} \chi(t) \approx \mu(t) + \sum_{k=1}^{K} \xi_k \, \phi_k(t) \end{equation}\]
La carta \(T^2\) se contruye sobre los scores de la siguiente manera:
\[\boldsymbol{\xi} = \left( \langle \chi - \mu, \phi_1 \rangle,\; \langle \chi - \mu, \phi_2 \rangle,\; \ldots,\; \langle \chi - \mu, \phi_K \rangle \right)^{'}\] Donde
\[\begin{align} \operatorname{Cov}(\xi_i, \xi_j) &= \operatorname{Cov}\big(\langle \chi - \mu, \phi_i \rangle,\; \langle \chi - \mu, \phi_j \rangle \big) \\ &= 0 \quad \text{si } i \neq j \\[6pt] \operatorname{Var}(\xi_i) &= \operatorname{Cov}(\xi_i, \xi_i) \\ &= \operatorname{Cov}\big(\langle \chi - \mu, \phi_i \rangle,\; \langle \chi - \mu, \phi_i \rangle \big) \\ &= \lambda_i \end{align}\]
Por lo tanto \(\Sigma_{\boldsymbol{\xi}} = \operatorname{diag}(\lambda_1, \lambda_2, \ldots, \lambda_K)\)
Sea \(x_1,x_2,..,x_n\) realizaciones de \(\chi\) tal que:
\[\boldsymbol{\xi}_j = \left( \langle x_j - \bar{x}, \phi_1 \rangle,\; \langle x_j - \bar{x}, \phi_2 \rangle,\; \ldots,\; \langle x_j - \bar{x}, \phi_K \rangle \right)^{'}.\]
Entonces el estadístico $ T_j^2$ se define como:
\[\begin{equation} T_j^2 = \boldsymbol{\xi}_j^{\top}\,\Sigma_{\boldsymbol{\xi}}^{-1}\,\boldsymbol{\xi}_j = \sum_{k=1}^{K} \frac{\xi_{jk}^2}{\lambda_k}, \end{equation}\]
El estadístico \(Q_j\) se define como:
\[ \begin{equation} Q_j = \| \varepsilon_j \|^2 =\| x_j -\widehat{x_j}\|^2= \left\| x_j - \bar{x} - \sum_{k=1}^{K} \xi_{jk}\,\phi_k \right\|^2, \end{equation} \]
En la Fase I se realiza el FPCA a partir de una muestra en control.
El nivel de significancia conjunto para las cartas de control \(T^2\) y \(Q\) se fija como \(\alpha\). Dado que se utilizan dos estadísticas de monitoreo, el nivel de significancia individual para cada carta se define como:
\[ \alpha' = 1 - (1 - \alpha)^{1/p}, \]
donde \(p = 2\) corresponde al número de cartas de control. De esta manera, se garantiza que el nivel de significancia global sea \(\alpha\). Nosotros fijaremos \(\alpha = 0.05\) como nivel de significancia conjunto.
Por otra parte, dado que la distribución exacta de los estadísticos \(T^2\) y \(Q\) no es conocida, los límites de control se estiman mediante un procedimiento de bootstrap. En particular, se considera el máximo de estos estadísticos en cada réplica bootstrap, y a partir de esta distribución empírica se calculan los cuantiles correspondientes de \(T^2\) y \(Q\), los cuales se utilizan para definir los limites.
\[ \text{LC}_{T^2} = q_{1-\alpha'}\left( {T_1^2}, \ldots, T_n^2 \right), \]
\[ \text{LC}_{Q} = q_{1-\alpha'}\left( Q_1, \ldots, Q_n \right), \]
Este enfoque permite evitar supuestos distribucionales y proporciona límites de control adaptados a la estructura de los datos.
Al final de la Fase I, se debe quedar con el conjunto de datos que represente un proceso en control.
Para una nueva observación \(x_{\text{new}}\), se calculan:
\[ T^2 = \sum_{k=1}^{K} \frac{\xi_k^2}{\lambda_k}, \qquad Q = \left\| x_{\text{new}} - \bar{x} - \sum_{k=1}^{K} \xi_k\,\phi_k \right\|^2. \]
Se considera fuera de control si:
\[ T^2 > \text{UCL}_{T^2} \quad \text{o} \quad Q > \text{UCL}_{Q}. \] Para la construcción de las cartas de control, se divide el conjunto de datos en dos subconjuntos: el \(80\%\) de las observaciones se utiliza en la Fase I, mientras que el \(20\%\) restante se reserva para la Fase II.
Una vez identificados los puntos de cambio, la construcción de las cartas de control se realizará mediante la partición del conjunto de datos, agrupando aquellas observaciones que comparten una misma estructura de media.
T2_Q.CL <- function(data, argvals, Lambda, Basis, ncomp, alpha, B, seed){
N <- dim(data)[2]
all_T2 <- vector("list", B)
all_Q <- vector("list", B)
fdParobj <- fdPar(fdobj = Basis, Lfdobj = 2, lambda = Lambda)
dt <- diff(argvals)
set.seed(seed)
for (b in seq_len(B)) {
# Remuestreo
idx <- sample(1:N, size = N, replace = TRUE)
samples <- data[, idx, drop = FALSE]
# Suavizado y PCA
smoothlist <- smooth.basis(argvals, samples, fdParobj)
fdObj_c <- center.fd(smoothlist$fd)
pcaObj <- pca.fd(fdObj_c, nharm = ncomp, centerfns = FALSE)
# Extraer scores y valores propios
Scores <- pcaObj$scores
Eigenval <- pcaObj$values[1:ncomp]
# T2
all_T2[[b]] <- rowSums((Scores %*% diag(1/Eigenval)) * Scores)
# Q
Xfd <- eval.fd(argvals, fdObj_c)
Phi <- eval.fd(argvals, pcaObj$harmonics)
Xhat <- Phi %*% t(Scores)
all_Q[[b]] <- apply((Xfd - Xhat)^2, 2, function(v){
sum((v[-1] + v[-length(v)]) / 2 * dt)
})
}
# Límites de control
CL_T2 <- quantile(unlist(all_T2), 1 - alpha)
CL_Q <- quantile(unlist(all_Q), 1 - alpha)
# Modelo con datos originales
smoothlist_f <- smooth.basis(argvals, data, fdParobj)
fdObj_c_f <- center.fd(smoothlist_f$fd)
pcaObj_f <- pca.fd(fdObj_c_f, nharm = ncomp, centerfns = FALSE)
Scores_f <- pcaObj_f$scores
Eigenval_f <- pcaObj_f$values[1:ncomp]
# T2 final
T2_final <- rowSums((Scores_f %*% diag(1/Eigenval_f)) * Scores_f)
# Q final
Xfd_f <- eval.fd(argvals, fdObj_c_f)
Phi_f <- eval.fd(argvals, pcaObj_f$harmonics)
Xhat_f <- Phi_f %*% t(Scores_f)
Q_final <- apply((Xfd_f - Xhat_f)^2, 2, function(v){
sum((v[-1] + v[-length(v)]) / 2 * dt)
})
return(list(
CL_T2 = as.numeric(CL_T2),
CL_Q = as.numeric(CL_Q),
T2 = T2_final,
Q = Q_final
))
}
clean_phaseI_T2Q <- function(data, p, argvals, Lambdas, Basis, ncomp = 2,
alpha = 0.025, B = 500, seed = 123,r,cp){
set.seed(seed)
data_variable <- data[ , , p][,cp]
n_total <- dim(data_variable)[2]
# Selección de Fase I
idx_phase1 <- 1:floor(r * n_total)
data_phase1 <- data_variable[, idx_phase1, drop = FALSE]
indices <- idx_phase1
removed_indices <- c()
Lambda <- Lambdas[p]
repeat {
res <- T2_Q.CL(
data = data_phase1,
argvals = argvals,
Lambda = Lambda,
Basis = Basis,
ncomp = ncomp,
alpha = alpha,
B = B,
seed = seed
)
# Identificar fuera de control
outliers <- which(res$T2 > res$CL_T2 | res$Q > res$CL_Q)
if(length(outliers) == 0) break
# Guardar cuáles estamos quitando
removed_indices <- c(removed_indices, indices[outliers])
# Actualizar datos y los índices que quedan
data_phase1 <- data_phase1[, -outliers, drop = FALSE]
indices <- indices[-outliers]
if(length(indices) <= ncomp) break
}
return(list(
cleaned_data = data_phase1,
final_result = res,
removed_indices = removed_indices,
n_removed = length(removed_indices),
indices = indices
))
}Se considerará una partición del conjunto de datos de la siguiente manera: la Fase I estará conformada por el grupo tres, correspondiente a las últimas observaciones que comparten una misma media. En este caso, la separación se realizará utilizando el 80% de los datos para la Fase I y el 20% restante para la Fase II. Los límites de control serán estimados mediante bootstrap.
Del último conjunto de datos, conformado por las últimas 84 observaciones que presentan la misma media, se seleccionó el 80%, correspondiente a 67 observaciones, para la Fase I. Dentro de este conjunto, 43 observaciones se utilizaron para la estimación de los límites de control, mientras que las 17 observaciones restantes se destinaron a la Fase II, con el propósito de realizar el monitoreo del proceso.
A partir de la Fase I, los límites de control obtenidos fueron 5.910113 para el estadístico \(T^2\) y 297.9885 para el estadístico \(Q\).
| Observación real | Valor T² | Valor Q | Fuera por T² | Fuera por Q |
|---|---|---|---|---|
| 254 | 1.311 | 450.092 | No | Sí |
| 255 | 0.502 | 424.033 | No | Sí |
| 259 | 2.871 | 317.485 | No | Sí |
| 266 | 0.644 | 475.134 | No | Sí |
| 268 | 26.473 | 1589.664 | Sí | Sí |
Sea \(\chi\) una función aleatoria definida en \(L^2\), asumimos \[ \mu_\chi(t) = E[\chi(t)], \] y \(\mathcal{C}_X\) el operador de covarianza \[ \mathcal{C}_\chi(\eta) = E\big[(X - \mu_X) \otimes (X - \mu_X)(\eta)\big]. \]
Sean \(\lambda_1, \lambda_2, \ldots\) los valores propios del operador y \(\psi_1, \psi_2, \ldots\) las funciones propias, que forman una base ortonormal de \(L^2\).
Usando la expansión de Karhunen–Loève de la función aleatoria \(X\): \[ \chi = \mu_\chi + \sum_{k=1}^{\infty} \xi_k \psi_k, \] donde \[ \xi_k = \langle \chi - \mu_X, \psi_k \rangle, \] los cuales son los scores de componentes principales, variables aleatorias no correlacionadas con media cero y varianza \(\lambda_k\).
El operador de covarianza se expresa como: \[ \mathcal{C}_\chi(x,y) = \sum_{j=1}^{\infty} \lambda_j \psi_j(x)\psi_j(y), \] usando el teorema de Mercer.
Para una función \(y\), el operador inverso se escribe como: \[ \mathcal{C}_\chi^{-1}(y) = \sum_{k=1}^{\infty} \frac{1}{\lambda_k} \langle y, \psi_k \rangle \psi_k. \]
Se define la distancia funcional de Mahalanobis entre \(\chi\) y \(\mu_\chi\) como: \[ d_{FM}^2 = \left\langle \mathcal{C}_\chi^{-1/2}(\chi - \mu_\chi), \mathcal{C}_\chi^{-1/2}(\chi - \mu_\chi) \right\rangle. \]
Expandiendo: \[ \mathcal{C}_\chi^{-1/2}(\chi - \mu_\chi) = \sum_{k=1}^{\infty} \frac{1}{\sqrt{\lambda_k}} \xi_k \psi_k, \] por lo que \[ d_{FM}^2 = \sum_{k=1}^{\infty} \frac{\xi_k^2}{\lambda_k}. \]
Equivalentemente, \[ d_{FM} = \left( \sum_{k=1}^{K} \frac{\xi_k^2}{\lambda_k} \right)^{1/2}. \]
Esta versión funcional de la distancia de Mahalanobis es una semi-distancia, ya que debido al truncamiento del operador de covarianza puede ocurrir que elementos distintos tengan distancia cero.
Para la construcción de las cartas de control se seguirá un procedimiento análogo al utilizado previamente en el contexto de FPCA, dividiendo el análisis en dos etapas: una fase de calibración (Fase I) y una fase de monitoreo (Fase II).
En la fase I se seleccionará un subconjunto de observaciones que se asume proviene de un proceso bajo control, es decir, datos que comparten la misma estructura de media. En este caso, se considerará nuevamente que el periodo estable corresponde a las observaciones finales, por lo que se tomará el 80% de los datos más recientes para conformar la Fase I, mientras que el 20% restante se destinará a la Fase II.
Esta partición permitirá estimar los límites de control a partir de un conjunto de referencia estable y, posteriormente, evaluar el comportamiento de las observaciones restantes. Manteniendo el mismo enfoque basado en bootstrap y en la distancia funcional de Mahalanobis.
DFM.CL <- function(data, argvals, Lambda, Basis, ncomp, alpha, B, seed){
N <- dim(data)[2]
all_DFM <- vector("list", B)
fdParobj <- fdPar(fdobj = Basis, Lfdobj = 2, lambda = Lambda)
set.seed(seed)
for (b in seq_len(B)) {
# Remuestreo bootstrap
idx <- sample(1:N, size = N, replace = TRUE)
samples <- data[, idx, drop = FALSE]
# Suavizado y FPCA
smoothlist <- smooth.basis(argvals, samples, fdParobj)
fdObj_c <- center.fd(smoothlist$fd)
pcaObj <- pca.fd(fdObj_c, nharm = ncomp, centerfns = FALSE)
# Scores y valores propios
Scores <- pcaObj$scores
Eigenval <- pcaObj$values[1:ncomp]
# Distancia funcional de Mahalanobis
all_DFM[[b]] <- sqrt(
rowSums(
(Scores^2) / matrix(
Eigenval,
nrow = nrow(Scores),
ncol = length(Eigenval),
byrow = TRUE
)
)
)
}
CL_DFM <- quantile(unlist(all_DFM), 1 - alpha)
# Modelo final con los datos originales
smoothlist_f <- smooth.basis(argvals, data, fdParobj)
fdObj_c_f <- center.fd(smoothlist_f$fd)
pcaObj_f <- pca.fd(fdObj_c_f, nharm = ncomp, centerfns = FALSE)
Scores_f <- pcaObj_f$scores
Eigenval_f <- pcaObj_f$values[1:ncomp]
# DFM final
DFM_final <- sqrt(
rowSums(
(Scores_f^2) / matrix(
Eigenval_f,
nrow = nrow(Scores_f),
ncol = length(Eigenval_f),
byrow = TRUE
)
)
)
return(list(
CL_DFM = as.numeric(CL_DFM),
DFM = DFM_final
))
}
clean_phaseI_DFM <- function(data, p, argvals,Lambdas,
Basis, ncomp = 2, alpha = 0.025, B = 500, seed = 123,
r, cp
) {
set.seed(seed)
data_variable <- data[, , p][, cp]
n_total <- dim(data_variable)[2]
# Selección aleatoria de Fase I
idx_phase1 <- 1:floor(r * n_total)
data_phase1 <- data_variable[, idx_phase1, drop = FALSE]
# Guardar índices originales
indices <- idx_phase1
removed_indices <- c()
Lambda <- Lambdas[p]
repeat {
# Estimar límite de control y DFM
res <- DFM.CL(
data = data_phase1,
argvals = argvals,
Lambda = Lambda,
Basis = Basis,
ncomp = ncomp,
alpha = alpha,
B = B,
seed = seed
)
# Identificar observaciones fuera de control
outliers <- which(res$DFM > res$CL_DFM)
# Si no hay outliers, detener
if (length(outliers) == 0) break
# Guardar índices eliminados
removed_indices <- c(
removed_indices,
indices[outliers]
)
# Eliminar observaciones fuera de control
data_phase1 <- data_phase1[, -outliers, drop = FALSE]
indices <- indices[-outliers]
# Evitar que queden muy pocas observaciones
if (length(indices) <= ncomp) break
}
return(list(
cleaned_data = data_phase1,
final_result = res,
removed_indices = removed_indices,
n_removed = length(removed_indices),
indices = indices
))
}
CL_DFM=clean_phaseI_DFM (data=wvl, p=2, argvals=t, Lambdas=Lambdas, Basis=Basis,B = 500,alpha=0.05,r=0.8,cp=185:268)Del último conjunto de datos, conformado por las últimas 84 observaciones que presentan la misma media, se seleccionó el 80%, correspondiente a 67 observaciones, para la Fase I. Dentro de este conjunto, 41 observaciones se utilizaron para la estimación de los límites de control, mientras que las 17 observaciones restantes se destinaron a la Fase II, con el propósito de realizar el monitoreo del proceso.
A partir de la Fase I, el límite de control obtenido fue 2.177862
| Observación real | Valor DFM |
|---|---|
| 267 | 2.357 |
| 268 | 5.145 |
Dado que el MFPCA se basa en la descomposición de Karhunen–Loève, se tiene que:
El vector funcional multivariante
\[\mathbf{X} = \big(X^{(1)}, \dots, X^{(p)}\big)\]
admite una representación finita de Karhunen–Loève si y solo si cada uno de los procesos univariantes
\(X^{(j)}, \quad j = 1, \dots, p\) admite una representación finita de Karhunen–Loève.
Esto implica que cada proceso univariante pertenece a un espacio de Hilbert, en particular al espacio de funciones aleatorias cuadrado integrables \(L^2(\mathcal{T})\), y que el operador de covarianza asociado es de tipo Hilbert–Schmidt, autoadjunto (simétrico) y definido positivo.
Para determinar las componentes principales funcionales multivariadas, se seguirá el enfoque propuesto por (Happ y Greven 2018), el cual permite realizar MFPCA incluso cuando los dominios de los procesos funcionales son diferentes.
En el presente caso, los dominios son iguales; no obstante, este enfoque resulta conveniente, ya que garantiza propiedades deseables del operador de covarianza del proceso funcional multivariado, tales como ser de tipo Hilbert–Schmidt, autoadjunto y definido positivo.
Xc1 <- center.fd(fit[[1]]$fd)
Xc2 <- center.fd(fit[[2]]$fd)
Xc3 <- center.fd(fit[[3]]$fd)
Xc4 <- center.fd(fit[[4]]$fd)
Xc5 <- center.fd(fit[[5]]$fd)
Xc6 <- center.fd(fit[[6]]$fd)
Xc7 <- center.fd(fit[[7]]$fd)Dadas muestras \(\mathbf{\chi_1}, \dots, \mathbf{\chi_N}\) de \(\mathbf{\chi}\), el procedimiento de estimación propuesto para MFPCA consta de cuatro pasos:
pcaObj1 <- pca.fd(Xc1, nharm = 2)
pcaObj2 <- pca.fd(Xc2, nharm = 2)
pcaObj3 <- pca.fd(Xc3, nharm = 2)
pcaObj4 <- pca.fd(Xc4, nharm = 2)
pcaObj5 <- pca.fd(Xc5, nharm = 2)
pcaObj6 <- pca.fd(Xc6, nharm = 2)
pcaObj7 <- pca.fd(Xc7, nharm = 2)pca_list <- list(pcaObj1, pcaObj2, pcaObj3, pcaObj4,
pcaObj5, pcaObj6, pcaObj7)
var_exp <- sapply(pca_list, function(obj) {
sum(obj$varprop)
})
tabla <- data.frame(
Proceso = paste0("$\\chi^{(", 1:7, ")}$"),
Varianza_Explicada = round(var_exp, 4)
)
knitr::kable(tabla, escape = FALSE)| Proceso | Varianza_Explicada |
|---|---|
| \(\chi^{(1)}\) | 0.9424 |
| \(\chi^{(2)}\) | 0.9633 |
| \(\chi^{(3)}\) | 0.9580 |
| \(\chi^{(4)}\) | 0.9645 |
| \(\chi^{(5)}\) | 0.9539 |
| \(\chi^{(6)}\) | 0.9416 |
| \(\chi^{(7)}\) | 0.9483 |
Definir la matriz \(\Xi \in \mathbb{R}^{N \times M_+}\), donde cada fila \(\Xi_i\) contiene todos los scores estimados para una observación: \(\Xi_i = \big(\hat{\xi}^{(1)}_{i,1}, \dots, \hat{\xi}^{(1)}_{i,M_1}, \dots, \hat{\xi}^{(p)}_{i,1}, \dots, \hat{\xi}^{(p)}_{i,M_p}\big).\) Un estimador \(\hat{Z} \in \mathbb{R}^{M_+ \times M_+}\) de la matriz de bloques \(Z\) se obtiene como \(\hat{Z} = (N - 1)^{-1} \Xi^\top \Xi.\)
Realizar un análisis espectral de matrices para \(\hat{Z}\), obteniendo valores propios \(\hat{\nu}_m\) y vectores propios ortonormales \(\hat{\mathbf{c}}_m\), \(m = 1, \dots, M_+\).
Los estimadores elemento a elemento de las funciones propias multivariantes asociadas a \(\hat{\nu}_m\) están dados por: \(\hat{\psi}^{(j)}_m(t_j) = \sum_{n=1}^{M_j} [\hat{\mathbf{c}}_m]^{(j)}_n \, \hat{\phi}^{(j)}_n(t_j), t_j \in \mathcal{T}_j,\) y los scores multivariantes pueden calcularse como:
Xc_list <- list(Xc1, Xc2, Xc3, Xc4, Xc5, Xc6, Xc7)
fit.funData <- list()
eigenfun.funData <- list()
scores_list <- list()
for (i in 1:7) {
fit.funData[[i]] <- fd2funData(Xc_list[[i]], argvals = grid)
eigenfun.funData[[i]] <- fd2funData(pca_list[[i]]$harmonics, argvals = grid)
scores_list[[i]] <- pca_list[[i]]$scores
}
mFData <- multiFunData(fit.funData)
uniExpansions <- list(list(type = "given",
functions = eigenfun.funData[[1]],
scores = scores_list[[1]], ortho = TRUE),
list(type = "given",
functions = eigenfun.funData[[2]],
scores = scores_list[[2]], ortho = TRUE),
list(type = "given",
functions = eigenfun.funData[[3]],
scores = scores_list[[3]], ortho = TRUE),
list(type = "given",
functions = eigenfun.funData[[4]],
scores = scores_list[[4]], ortho = TRUE),
list(type = "given",
functions = eigenfun.funData[[5]],
scores = scores_list[[5]], ortho = TRUE),
list(type = "given",
functions = eigenfun.funData[[6]],
scores = scores_list[[6]], ortho = TRUE),
list(type = "given",
functions = eigenfun.funData[[7]],
scores = scores_list[[7]], ortho = TRUE))
MFPCAObj <- MFPCA(mFData,
M = 3,
uniExpansions = uniExpansions,
fit = TRUE,
verbose = TRUE)Calculating univariate basis expansions (02:28:56)
Calculating MFPCA (02:28:56)par(mfrow = c(2,4), mar = c(3,3,2,1))
cols <- 1:3
for (j in 1:7) {
matplot(MFPCAObj$functions[[j]]@argvals[[1]],
t(MFPCAObj$functions[[j]]@X),
type = "l", lty = 1,
col = cols,
xlab = "t", ylab = "",
main = paste("Proceso", j))
legend("topright",
legend = c(expression(psi[1]),
expression(psi[2]),
expression(psi[3])),
col = cols, lty = 1, lwd = 2,
cex = 0.8, bty = "n")
}
Xc_list <- list(Xc1, Xc2, Xc3, Xc4, Xc5, Xc6, Xc7)
par(mfrow = c(2,4), mar = c(3,3,2,1)) # 7 plots (2x4)
for (j in 1:7) {
plot(Xc_list[[j]],
xlab = "t", ylab = "",
main = paste("Proceso", j))
}
\(\Sigma_{\chi^{(1)}\chi^{(1)}}(t, s) = (N - 1)^{-1} \sum_{i=1}^{N} \{x_{i}^{(1)}(t) - \bar{x}^{(1)}(t)\}\{x_{i}^{(1)}(s) - \bar{x}^{(1)}(s)\}\)
wvl_c1 <- sweep(wvl[,,1], 1, rowMeans(wvl[,,1]), "-")
wvl_c2 <- sweep(wvl[,,2], 1, rowMeans(wvl[,,2]), "-")
sigma11 = (wvl_c1 %*% t(wvl_c1))/(268-1)
sigma12 = (wvl_c1 %*% t(wvl_c2))/(268-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 = ~t, y = ~t, z = ~sigma11) %>%
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 Σ11"),
scene = list(
xaxis = list(title = "Longitud de onda (nm)"),
yaxis = list(title = "Longitud de onda (nm)"),
zaxis = list(
title = "Valor",
range = c(0,4567),
tickformat = ".0f"
)
)
)
p\(\Sigma_{\chi^{(1)}\chi^{(2)}}(t, s) = (N - 1)^{-1} \sum_{i=1}^{N} \{x_{i}^{(1)}(t) - \bar{x}^{(1)}(t)\}\{x_{i}^{(2)}(s) - \bar{x}^{(2)}(s)\}\)
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 = ~t, y = ~t, z = ~sigma11) %>%
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 Σ12"),
scene = list(
xaxis = list(title = "Longitud de onda (nm)"),
yaxis = list(title = "Longitud de onda (nm)"),
zaxis = list(
title = "Valor",
range = c(0,827),
tickformat = ".0f"
)
)
)
pvar_fd <- var.fd(fit[[1]]$fd)
Cmat <- eval.bifd(t, t, var_fd)
plot(t,diag(Cmat))
var_exp <- MFPCAObj14[["values"]] / sum(MFPCAObj14[["values"]])
cum_var_exp <- cumsum(var_exp)plot(cum_var_exp[1:8], type = "b",
xlab = "Número de componentes",
ylab = "Varianza explicada acumulada",
ylim = c(0, 1.2),
pch = 19)
abline(h = 0.9, col = "red", lty = 2)
text(1:length(cum_var_exp),
cum_var_exp,
labels = round(cum_var_exp, 3),
pos = 3,
cex = 0.8)
load("MFPCAObj2.RData")par(mfrow = c(2,4), mar = c(3,3,2,1))
cols <- 1:2
for (j in 1:7) {
matplot(MFPCAObj2$functions[[j]]@argvals[[1]],
t(MFPCAObj2$functions[[j]]@X),
type = "l", lty = 1,
col = cols,
xlab = "t", ylab = "",
main = paste("Proceso", j))
legend("topright",
legend = c(expression(psi[1]),
expression(psi[2])),
col = cols, lty = 1, lwd = 2,
cex = 0.8, bty = "n")
}
values <- MFPCAObj14[["values"]]
plot(log(values), type = "b",
pch = 19, lwd = 2,
col = "darkred",
ylim = c(0,14),
xlab = "Número de componentes",
ylab = "Valor propio (escala log)")
grid()
round(colMeans(pcaObj1$scores),4)[1] 0 0round(colMeans(pcaObj2$scores),4)[1] 0 0cor_matrix <- function(obj1, obj2){
matrix(c(
cor(obj1$scores[,1], obj2$scores[,1]),
cor(obj1$scores[,1], obj2$scores[,2]),
cor(obj1$scores[,2], obj2$scores[,1]),
cor(obj1$scores[,2], obj2$scores[,2])
), nrow = 2, byrow = TRUE)
}
cor_scores_12 <- cor_matrix(pcaObj1,pcaObj2)
rownames(cor_scores_12) <- c("E11", "E21")
colnames(cor_scores_12) <- c("E12", "E22")
knitr::kable(cor_scores_12, digits = 3)| E12 | E22 | |
|---|---|---|
| E11 | 0.800 | 0.405 |
| E21 | 0.566 | -0.329 |
pca_list <- list(pcaObj1, pcaObj2, pcaObj3, pcaObj4,
pcaObj5, pcaObj6, pcaObj7)
n <- length(pca_list)
cor_list <- list()
for(i in 1:n){
for(j in i:n){
cor_list[[paste0("Proc",i,"-Proc",j)]] <- cor_matrix(pca_list[[i]], pca_list[[j]])
}
}# inicializar matriz grande (2 PCs por proceso)
cor_big <- matrix(NA, nrow = 2*n, ncol = 2*n)
# nombres
labels <- unlist(lapply(1:n, function(i) paste0("P", i, "_PC", 1:2)))
rownames(cor_big) <- labels
colnames(cor_big) <- labelsfor(i in 1:n){
for(j in 1:n){
# matriz 2x2 entre proceso i y j
block <- cor_matrix(pca_list[[i]], pca_list[[j]])
# posiciones en la matriz grande
rows <- ((i-1)*2 + 1):(i*2)
cols <- ((j-1)*2 + 1):(j*2)
cor_big[rows, cols] <- round(block,2)
}
}cor_long <- melt(cor_big)
ggplot(cor_long, aes(Var1, Var2, fill = value)) +
geom_tile() +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limits = c(-1,1)) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(x = "", y = "", fill = "Correlación")
cor_scores_multi <- cor(MFPCAObj14$scores)
cor_long <- melt(cor_scores_multi)
ggplot(cor_long, aes(Var1, Var2, fill = value)) +
geom_tile() +
geom_text(aes(label = round(value, 2)), size = 3) +
scale_fill_gradient2(low = "blue", high = "red", mid = "white",
midpoint = 0, limits = c(-1,1)) +
theme_minimal() +
labs(x = "Componentes", y = "Componentes", fill = "Correlación") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
df <- data.frame(
PC1 = MFPCAObj2$scores[,1],
PC2 = MFPCAObj2$scores[,2]
)
ggplot(df, aes(x = PC1, y = PC2)) +
geom_point(alpha = 0.6, size = 2, color = "steelblue") +
theme_minimal() +
labs(
title = "Dispersión de los scores (PC1 vs PC2)",
x = "Componente principal 1",
y = "Componente principal 2"
) +
geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
geom_vline(xintercept = 0, linetype = "dashed", color = "gray")
En el caso multivariado, se tiene la descomposición de Karhunen–Loève, de forma que para el vector de curvas aleatorias, $ L^2( _1 ) L^2( _p) $, se tiene que
\[ \mathbf{X}(\mathbf{t}) = \boldsymbol{\mu}(\mathbf{t}) + \sum_{k=1}^\infty \xi_k \boldsymbol{\psi}_k (\mathbf{t}) \] con \(\mathbf{t} = (t_1,\ldots,t_p)^t\). En este caso, se tiene un vector de funciones de tamaño \(p=7\), donde el dominio es el mismo para cada proceso, \(\mathcal{T_1}=\cdots=\mathcal{T_7} = [275, 560]\). De esta forma, se aproxima cada una de las observaciones muestrales de la siguiente manera:
\[ \mathbf{X}_i(\mathbf{t}) \approx \hat {\boldsymbol{\mu}}(\mathbf{t}) + \sum_{k=1}^K \hat{\xi}_{ik} \hat{\boldsymbol{\psi}}_k (\mathbf{t}) \]
X_list <- list(Xc1, Xc2, Xc3, Xc4, Xc5, Xc6, Xc7)
n_proc <- length(X_list)
for(j in 1:n_proc){
# aproximaciones
aprox3 <- MFPCAObj2[["fit"]][[j]]@X[1:3,]
# dominio
t <- MFPCAObj2[["functions"]][[j]]@argvals[[1]]
# funciones originales
muestra3 <- t(eval.fd(t, X_list[[j]])[,1:3])
# gráfico
par(mfrow = c(1,3))
for(i in 1:3){
plot(t, muestra3[i,], type = "l", lwd = 2,
col = "blue",
ylim = range(c(muestra3, aprox3)),
xlab = "t", ylab = "Valor",
main = paste("Proceso", j, "- Curva", i))
lines(t, aprox3[i,], col = "red", lwd = 2, lty = 2)
legend("topright",
legend = c("Original", "Aprox"),
col = c("blue", "red"),
lty = c(1,2),
bty = "n")
}
}
En la sección 4 del artículo de (Happ y Greven 2018) se presentan varias simulaciones para medir el desempeño de su propuesta. Se tienen tres escenarios:
En cada caso, se tienen 100 datasets con \(N=250\) observaciones: \[ \mathbf{x}_i(\mathbf{t}) = \sum_{m=1}^{M} \rho_{im}\psi_m(\mathbf{t}) + \boldsymbol{\varepsilon}_i(\mathbf{t}), \qquad \boldsymbol{\varepsilon}_i(\mathbf{t}) \overset{\mathrm{iid}}{\sim} \mathrm{N}_p(0, \sigma^2 \mathbf{I}) \] con \(\mathbf{t} \in \mathcal{T}\) e \(i=1,\ldots, N = 250\).
Además se tienen otras configuraciones:
\[ \nu_m^{\text{exp}} = \exp{(-(m+1)/2)}, \qquad \nu_m^{\text{lin}} = (M+1-m)/M \]
La precisión de las estimaciones se mide con el error relativo
\[ \text{Err}(\hat{\nu}_m) = \frac{(\nu_m - \hat{\nu}_m)^2}{\nu_m^2} \]
donde \(\nu_m\) es el valor propio verdadero y \(\hat{\nu}_m\) es el valor propio estimado. Valores pequeños de esta medida indican una buena estimación del valor propio correspondiente. Por otro lado, la precisión en la estimación de las autofunciones se evalúa mediante \[ \text{Err}(\hat{\psi}_m) = ||| \psi_m - \hat{\psi}_m |||^2, \]
donde \(\psi_m\) corresponde a la función propia verdadera y \(\hat{\psi}_m\) a su estimación. Debido a que las funciones propias están definidas salvo un cambio de signo, se ajusta la orientación de la función propia estimada cuando \[ \langle\langle \psi_m, \hat{\psi}_m \rangle\rangle < 0. \]
En ese caso, se reemplaza \(\hat{\psi}_m\) por \(-\hat{\psi}_m\), evitando considerar como error una diferencia que solo se debe a la orientación de la componente. Además de evaluar autovalores y autofunciones, las autoras analizan la calidad de reconstrucción de las observaciones originales. Para ello utilizan el error medio relativo cuadrático, definido como
\[ \text{MRSE} = \frac{1}{N} \sum_{i=1}^{N} \frac{|||x_i - \hat{x}_i|||^2}{|||x_i|||^2}. \]
La reconstrucción de cada observación se obtiene mediante
\[ \hat{x}_i = \sum_{m=1}^{M} \hat{\rho}_{i,m}\hat{\psi}_m. \]
En la primera configuración se tienen en cuenta las primeras \(M=8\) funciones de la base de Fourier en el intervalo \([0,1]\) para ambos procesos funcionales que forman el proceso bivariado de forma que \(\mathcal{T_1}=\mathcal{T_2}=[0,1]\).
Las obsevaciones se muestrean de un grid equidistante de \(S_1=S_2=100\) puntos y se usan \(M_1=M_2=8\) componentes principales funcionales univariadas para estimar las \(M=8\) componentes principales funcionales multivariadas. El método propuesto por los autores se compara con el método implementado en fda, que se denotará como \(\mathrm{MFPCA_{RS}}\), donde se realiza un suavizamiento previo con \(K=15\) splines cúbicos. La función pca.fd(), a diferencia de MFPCA(), devuelve la estimación de dos scores por observación bivariada y función propia estimada, de forma que se suman para poder compararlos con la nueva metodología.
La simulación se realiza mediante la función simMultiFunData() del paquete funData.
Sim_Setting1 <- function(eValType = c("exponential", "linear"),
s2 = 0,
seed = NULL,
N = 250,
M = 8,
argvals = argvalsSetting1,
argvalsfd,
uniExpansions = uniExpansions,
B){
set.seed(seed)
## MFPCA
MRSE <- numeric(B)
ErrValues <- matrix(NA, nrow = B, ncol = M)
ErrFuns <- matrix(NA, nrow = B, ncol = M)
## fda
MRSE_RS <- numeric(B)
ErrValuesRS <- matrix(NA, nrow = B, ncol = M)
ErrFunsRS <- matrix(NA, nrow = B, ncol = M)
for (i in 1:B) {
SimData <- simMultiFunData(type = "split",
argvals = argvalsSetting1,
M = M,
eFunType = "Fourier",
eValType = eValType,
N = N)
TrueValues <- SimData$trueVals
TrueFunctions <- SimData$trueFuns
if (s2 > 0) {
SimData2 <- addError(SimData$simData, sd = sqrt(s2))
} else {
SimData2 <- SimData$simData
}
MFPCAObj <- MFPCA(SimData2,
M = M,
uniExpansions = uniExpansions,
fit = TRUE)
FittedValues <- MFPCAObj$values
FittedFunctions <- MFPCAObj$functions
fits <- MFPCAObj$fit
# MRSE
num <- norm(SimData2 - fits, squared = TRUE)
den <- norm(SimData2, squared = TRUE)
MRSE[i] <- mean(num / den)
# Eigenvalue relative error
ErrValues[i, ] <- (TrueValues - FittedValues) ^ 2 / (TrueValues ^ 2)
# Eigenfunction relative error
ErrFunsA <- norm(TrueFunctions - FittedFunctions, squared = TRUE)
ErrFunsB <- norm(TrueFunctions + FittedFunctions, squared = TRUE)
ErrFunsC <- norm(TrueFunctions, squared = TRUE)
ErrFuns[i, ] <- pmin(ErrFunsA, ErrFunsB)
### fda package
fdData1 <- aperm(SimData2@.Data[[1]]@X, c(2, 1))
fdData2 <- aperm(SimData2@.Data[[2]]@X, c(2, 1))
fdData <- array(NA, dim = c(nrow(fdData1), ncol(fdData1), 2))
fdData[ , , 1] <- fdData1
fdData[ , , 2] <- fdData2
fdarange <- range(argvalsfd)
fdabasis = create.bspline.basis(rangeval = fdarange,
nbasis = 15,
norder = 4)
fdatime <- argvalsfd
fdafd <- smooth.basis(fdatime, fdData, fdabasis)$fd
nharm <- M
fdapcaList <- pca.fd(fdafd, nharm)
fdaValues <- fdapcaList$values[1:M]
harmEval <- eval.fd(argvalsfd, fdapcaList$harmonics)
fdaFuns <- multiFunData(list(funData(argvals = list(argvalsfd),
X = t(harmEval[, , 1])),
funData(argvals = list(argvalsfd),
X = t(harmEval[, , 2]))))
sgn <- ifelse(scalarProduct(TrueFunctions, fdaFuns) < 0, -1, 1)
for (j in 1:2) {
fdaFuns[[j]]@X <- fdaFuns[[j]]@X * sgn
}
### MRSE
meanEval <- eval.fd(argvalsfd, fdapcaList$meanfd)
scores_RS <- apply(fdapcaList$scores[, , , drop = FALSE], c(1, 2), sum)
fits_fd <- array(0, dim = c(length(argvalsfd), N, 2))
for (j in 1:2) {
for (m in 1:M) {
fits_fd[, , j] <- fits_fd[, , j] + outer(harmEval[, m, j], scores_RS[, m])
}
fits_fd[, , j] <- sweep(fits_fd[, , j], 1, meanEval[, , j], "+")
}
fits_RS <- multiFunData(list(funData(argvals = list(argvalsfd),
X = t(fits_fd[, , 1])),
funData(argvals = list(argvalsfd),
X = t(fits_fd[, , 2]))))
numRS <- norm(SimData2 - fits_RS, squared = TRUE)
denRS <- norm(SimData2, squared = TRUE)
MRSE_RS[i] <- mean(numRS / denRS)
# Eigenvalue relative error
ErrValuesRS[i, ] <- (TrueValues - fdaValues) ^ 2 / (TrueValues ^ 2)
# Eigenfunction relative error
ErrFunsARS <- norm(TrueFunctions - fdaFuns, squared = TRUE)
ErrFunsBRS <- norm(TrueFunctions + fdaFuns, squared = TRUE)
ErrFunsCRS <- norm(TrueFunctions, squared = TRUE)
ErrFunsRS[i, ] <- pmin(ErrFunsARS, ErrFunsBRS)
}
return(list(MRSE = MRSE,
ErrValues = ErrValues,
ErrFuns = ErrFuns,
MRSERS = MRSE_RS,
ErrValuesRS = ErrValuesRS,
ErrFunsRS = ErrFunsRS))
}A continuación, se define el número de observaciones a simular, el error, el dominio sobre el cual se generan las observaciones y el número de componentes principales univariadas a estimar. Además, se muestra un ejemplo de un dataset generado con valores propios que decaen exponencialmente en sin y con error de medición en la Figura 37 y Figura 38. De igual forma, en la Figura 39 se muestran las funciones propias simuladas.
N <- 250
s2 <- 0.25 # measurement error
argvalsSetting1 <- list(seq(0, 1, length = 100), # 100 sampling points
seq(0, 1, length = 100))
argvalsfd <- seq(0, 1, length = 100)
uniExpansions <- list(list(type = "uFPCA", npc = 8),
list(type = "uFPCA", npc = 8))
M <- 8
SimDataS1 <- simMultiFunData(type = "split",
argvals = argvalsSetting1,
M = M,
eFunType = "Fourier",
eValType = "exponential",
N = N)
SimDataS1s2 <- addError(SimDataS1$simData, sd = sqrt(s2))
A continuación, se realiza la simulación de esta primera configuración \(B=100\) veces con ambos tipos de decaimiento de los valores propios (exponencial y lineal) y errores de medida (\(\sigma^2=0\) y \(\sigma^2=0.25\)).
S1_exp_s0 <- Sim_Setting1(eValType = "exponential",
s2 = 0,
seed = 1,
N = N,
M = 8,
argvals = argvalsSetting1,
argvalsfd = argvalsfd,
uniExpansions = uniExpansions,
B = 100)
S1_exp_s <- Sim_Setting1(eValType = "exponential",
s2 = s2,
seed = 1,
N = N,
M = 8,
argvals = argvalsSetting1,
argvalsfd = argvalsfd,
uniExpansions = uniExpansions,
B = 100)
S1_lin_s0 <- Sim_Setting1(eValType = "linear",
s2 = 0,
seed = 1,
N = N,
M = 8,
argvals = argvalsSetting1,
argvalsfd = argvalsfd,
uniExpansions = uniExpansions,
B = 100)
S1_lin_s <- Sim_Setting1(eValType = "linear",
s2 = s2,
seed = 1,
N = N,
M = 8,
argvals = argvalsSetting1,
uniExpansions = uniExpansions,
argvalsfd = argvalsfd,
B = 100)Los resultados se presentan en la Figura 40 y la Tabla 4.
| Configuración | Método | \(\sigma^2 = 0\), \(\nu_m^{exp}\) | \(\sigma^2 = 0\), \(\nu_m^{lin}\) | \(\sigma^2 = 0.25\), \(\nu_m^{exp}\) | \(\sigma^2 = 0.25\), \(\nu_m^{lin}\) |
|---|---|---|---|---|---|
| 1 | MFPCA | 1e-04 | 1e-04 | 0.2165 | 0.1246 |
| 1 | MFPCARS | 0e+00 | 0e+00 | 0.2150 | 0.1238 |
Como se puede ver, la propuesta presentada en el paper compite bastante con método de Ramsay y Silverman que se encuentra implementado en el paquete fda, dando resultados bastante similares. En el caso de los valores propios con decaimiento exponencial, en los casos en los que no hay error de medición se tienen los resultados más parejos, mientras que cuando se incorpora el error de medición, las estimaciones de los últimios valores propios tienen un mayor error relativo. En el caso de las funciones propias, las estimaciones tienen errores muy similares.
Por otro lado, el MRSE nos permite verificar que el error de las reconstrucciones mediante la descomposición de Karhúnen–Loéve truncada es similar, con diferencia hasta el tercer decimal.
Para la segunda configuración, se evalúa el rendimiento de la nueva propuesta para datos dispersos con las primeras \(M=8\) funciones en la base de Fourier. Se realiza la comparación sin dispersión (Full data), dispersión media (Medium sparsity) y dispersión alta (High sparsity). A diferencia de la configuración anterior, se consideran realizaciones de un proceso funcional trivariado con distintos dominios (\(\mathcal{T}_1=[-1,-0.5]\), \(\mathcal{T}_2=[0,1]\), \(\mathcal{T}_3=[1.5,2]\)) y número de puntos discretizados distintos (\(S_1=S3=50\) y \(S_2=100\)). Además, para cada opción de sparsity se tienen las siguientes aclaraciones:
El MFPCA se calcula dependiendo de la dispersión. Para las simulaciones de datos completos y dispersión media se usan \(M_1=M_2=M_3=8\) componentes principales funcionales univariadas. Con dispersión alta, se usan \(M_1=M_3=3\) y \(M_2=5\) por estabilidad.
A continuación, en la Figura 41 y la Figura 42 se muestra un ejemplo de un dataset simulado con un decaimiento exponencial de los valores propios. Al igual que en el caso anterior, se hace uso de la función simMultiFunData(). Además, en la Figura 43 se muestran las funciones propias simuladas.
N <- 250
s2 <- 0.25 # measurement error
T1 <- seq(-1, -0.5, length = 50)
T2 <- seq(0, 1, length = 100)
T3 <- seq(1.5, 2, length = 50)
M <- 8
Ms <- c(8, 8, 8)
MsHigh <- c(3, 5, 3)
argvalsSetting2 <- list(T1, T2, T3)
M1 <- Ms[1]
M2 <- Ms[2]
M3 <- Ms[3]
M1High <- MsHigh[1]
M2High <- MsHigh[2]
M3High <- MsHigh[3]
uniExpansions <- list(list(type = "uFPCA", npc = M1),
list(type = "uFPCA", npc = M2),
list(type = "uFPCA", npc = M3))
uniExpansionsHigh <- list(list(type = "uFPCA", npc = M1High),
list(type = "uFPCA", npc = M2High),
list(type = "uFPCA", npc = M3High))
SimDataS2 <- simMultiFunData(type = "split",
argvals = argvalsSetting2,
M = M,
eFunType = "Fourier",
eValType = "exponential",
N = N)
SimDataS2s2 <- addError(SimDataS2$simData, sd = sqrt(s2))
Tlens <- unlist(lapply(argvalsSetting2, length))
minObs <- round(0.3 * Tlens)
maxObs <- round(0.5 * Tlens)
SimDataS2s2Med <- sparsify(SimDataS2$simData, minObs = minObs, maxObs = maxObs)
Para la simulación, primero se definen unas funciones para cada tipo de dispersión a considerar en la simulación.
FullSparsity <- function(Tlens, SimData, SimData2, M, uniExpansions){
TrueValues <- SimData$trueVals
TrueFunctions <- SimData$trueFuns
TrueObs <- SimData2
SimData3 <- SimData2
MFPCAObj <- MFPCA(SimData3,
M = M,
uniExpansions = uniExpansions,
fit = TRUE)
FittedValues <- MFPCAObj$values
FittedFunctions <- MFPCAObj$functions
fits <- MFPCAObj$fit
# MRSE
num <- norm(TrueObs - fits, squared = TRUE)
den <- norm(TrueObs, squared = TRUE)
MRSE <- mean(num / den)
# Eigenvalue relative error
ErrValues <- (TrueValues - FittedValues) ^ 2 / (TrueValues ^ 2)
# Eigenfunction relative error
ErrFunsA <- norm(TrueFunctions - FittedFunctions, squared = TRUE)
ErrFunsB <- norm(TrueFunctions + FittedFunctions, squared = TRUE)
ErrFunsC <- norm(TrueFunctions, squared = TRUE)
ErrFuns <- pmin(ErrFunsA, ErrFunsB)
return(list(MRSE = MRSE,
ErrValues = ErrValues,
ErrFuns = ErrFuns))
}
MediumSparsity <- function(Tlens, SimData, SimData2, M, uniExpansions){
TrueObs <- SimData2
TrueValues <- SimData$trueVals
TrueFunctions <- SimData$trueFuns
minObs <- round(0.3 * Tlens)
maxObs <- round(0.5 * Tlens)
SimData3 <- sparsify(SimData2, minObs = minObs, maxObs = maxObs)
MFPCAObj <- MFPCA(SimData3,
M = M,
uniExpansions = uniExpansions,
fit = TRUE)
FittedValues <- MFPCAObj$values
FittedFunctions <- MFPCAObj$functions
fits <- MFPCAObj$fit
# MRSE
num <- norm(TrueObs - fits, squared = TRUE)
den <- norm(TrueObs, squared = TRUE)
MRSE <- mean(num / den)
# Eigenvalue relative error
ErrValues <- (TrueValues - FittedValues) ^ 2 / (TrueValues ^ 2)
# Eigenfunction relative error
ErrFunsA <- norm(TrueFunctions - FittedFunctions, squared = TRUE)
ErrFunsB <- norm(TrueFunctions + FittedFunctions, squared = TRUE)
ErrFunsC <- norm(TrueFunctions, squared = TRUE)
ErrFuns <- pmin(ErrFunsA, ErrFunsB)
return(list(MRSE = MRSE,
ErrValues = ErrValues,
ErrFuns = ErrFuns))
}
HighSparsity <- function(Tlens, SimData, SimData2, M, uniExpansions, MsHigh){
TrueObs <- SimData2
TrueValues <- SimData$trueVals
TrueFunctions <- SimData$trueFuns
minObs <- pmax(round(0.1 * Tlens), MsHigh + 1)
maxObs <- round(0.3 * Tlens)
SimData3 <- sparsify(SimData2, minObs = minObs, maxObs = maxObs)
MFPCAObj <- MFPCA(SimData3,
M = M,
uniExpansions = uniExpansions,
fit = TRUE)
FittedValues <- MFPCAObj$values
FittedFunctions <- MFPCAObj$functions
fits <- MFPCAObj$fit
# MRSE
num <- norm(TrueObs - fits, squared = TRUE)
den <- norm(TrueObs, squared = TRUE)
MRSE <- mean(num / den)
# Eigenvalue relative error
ErrValues <- (TrueValues - FittedValues) ^ 2 / (TrueValues ^ 2)
# Eigenfunction relative error
ErrFunsA <- norm(TrueFunctions - FittedFunctions, squared = TRUE)
ErrFunsB <- norm(TrueFunctions + FittedFunctions, squared = TRUE)
ErrFunsC <- norm(TrueFunctions, squared = TRUE)
ErrFuns <- pmin(ErrFunsA, ErrFunsB)
return(list(MRSE = MRSE,
ErrValues = ErrValues,
ErrFuns = ErrFuns))
}La función para realizar las simulaciones se presenta a continuación.
Sim_Setting2 <- function(eValType = c("exponential", "linear"),
s2 = 0,
seed = NULL,
N = 250,
M,
Ms,
MsHigh,
argvals = argvalsSetting2,
B){
M1 <- Ms[1]
M2 <- Ms[2]
M3 <- Ms[3]
M1High <- MsHigh[1]
M2High <- MsHigh[2]
M3High <- MsHigh[3]
uniExpansions <- list(list(type = "uFPCA", npc = M1),
list(type = "uFPCA", npc = M2),
list(type = "uFPCA", npc = M3))
uniExpansionsHigh <- list(list(type = "uFPCA", npc = M1High),
list(type = "uFPCA", npc = M2High),
list(type = "uFPCA", npc = M3High))
Tlens <- unlist(lapply(argvals, length))
MRSE_Full <- numeric(B)
ErrValues_Full <- matrix(NA, nrow = B, ncol = M)
ErrFuns_Full <- matrix(NA, nrow = B, ncol = M)
MRSE_Medium <- numeric(B)
ErrValues_Medium <- matrix(NA, nrow = B, ncol = M)
ErrFuns_Medium <- matrix(NA, nrow = B, ncol = M)
MRSE_High <- numeric(B)
ErrValues_High <- matrix(NA, nrow = B, ncol = M)
ErrFuns_High <- matrix(NA, nrow = B, ncol = M)
set.seed(seed)
## MFPCA
for (i in 1:B) {
SimData <- simMultiFunData(type = "split",
argvals = argvals,
M = M,
eFunType = "Fourier",
eValType = eValType,
N = N)
if (s2 > 0) {
SimData2 <- addError(SimData$simData, sd = sqrt(s2))
} else {
SimData2 <- SimData$simData
}
Full <- FullSparsity(Tlens = Tlens,
SimData = SimData,
SimData2 = SimData2,
M = M,
uniExpansions = uniExpansions)
Medium <- MediumSparsity(Tlens = Tlens,
SimData = SimData,
SimData2 = SimData2,
M = M,
uniExpansions = uniExpansions)
High <- HighSparsity(Tlens = Tlens,
SimData = SimData,
SimData2 = SimData2,
M = M,
uniExpansions = uniExpansionsHigh,
MsHigh = MsHigh)
MRSE_Full[i] <- Full$MRSE
ErrValues_Full[i, ] <- Full$ErrValues
ErrFuns_Full[i, ] <- Full$ErrFuns
MRSE_Medium[i] <- Medium$MRSE
ErrValues_Medium[i, ] <- Medium$ErrValues
ErrFuns_Medium[i, ] <- Medium$ErrFuns
MRSE_High[i] <- High$MRSE
ErrValues_High[i, ] <- High$ErrValues
ErrFuns_High[i, ] <- High$ErrFuns
}
MRSE = list(Full = MRSE_Full,
Medium = MRSE_Medium,
High = MRSE_High)
ErrValues = list(Full = ErrValues_Full,
Medium = ErrValues_Medium,
High = ErrValues_High)
ErrFuns <- list(Full = ErrFuns_Full,
Medium = ErrFuns_Medium,
High = ErrFuns_High)
return(list(MRSE = MRSE,
ErrValues = ErrValues,
ErrFuns = ErrFuns))
}Y se realizan las respectivas simulaciones \(B=100\) veces, con ambos tipos de decaimiento de los valores propios, errores de medición y también los distintos niveles de sparsity.
N <- 250
s2 <- 0.25 # measurement error
T1 <- seq(-1, -0.5, length = 50)
T2 <- seq(0, 1, length = 100)
T3 <- seq(1.5, 2, length = 50)
M <- 8
Ms <- c(8, 8, 8)
MsHigh <- c(3, 5, 3)
argvalsSetting2 <- list(T1, T2, T3)
S2_exp_s0 <- Sim_Setting2(eValType = "exponential",
s2 = 0,
seed = 2,
N = N,
M = M,
Ms = Ms,
MsHigh = MsHigh,
argvals = argvalsSetting2,
B = 100)
S2_exp_s <- Sim_Setting2(eValType = "exponential",
s2 = s2,
seed = 2,
N = 250,
M = 8,
Ms = rep(8, 3),
MsHigh = c(3, 5, 3),
argvals = argvalsSetting2,
B = 100)
S2_lin_s0 <- Sim_Setting2(eValType = "linear",
s2 = 0,
seed = 2,
N = 250,
M = 8,
Ms = rep(8, 3),
MsHigh = c(3, 5, 3),
argvals = argvalsSetting2,
B = 100)
S2_lin_s <- Sim_Setting2(eValType = "linear",
s2 = s2,
seed = 2,
N = 250,
M = 8,
Ms = rep(8, 3),
MsHigh = c(3, 5, 3),
argvals = argvalsSetting2,
B = 100)Los resultados se presentan en la Figura 44.
| Método | \(\sigma^2 = 0\), \(\nu_m^{exp}\) | \(\sigma^2 = 0\), \(\nu_m^{lin}\) | \(\sigma^2 = 0.25\), \(\nu_m^{exp}\) | \(\sigma^2 = 0.25\), \(\nu_m^{lin}\) |
|---|---|---|---|---|
| Full data | 0.0062 | 0.0083 | 21.4450 | 12.3222 |
| Medium sparsity | 0.3253 | 0.2987 | 23.3041 | 13.5365 |
| High sparsity | 2.5808 | 2.0033 | 27.9791 | 16.9083 |
Como es posible ver en la Figura 44, para el escenario con los valores propios decayendo exponencialmente, hay mayor dificultar para estimar tanto los valores y funciones propias en el caso en el que hay una alta dispersión, teniendo los mayores errores relativos. Esto también ocurre en el caso de los valores propios con decaimiento lineal, pero se presentan menos dificultades en el proceso de estimación. Se puede ver que tanto en el escenario de los datos completos como de dipsersión media, se obtienen resultados muy similares. En general, el método mantiene un buen desempeño con datos completos y con dispersión moderada. Cuando la dispersión es alta, el error de reconstrucción aumenta, pero las primeras componentes principales siguen siendo estimadas de manera razonable. Esto es relevante porque en aplicaciones reales los datos funcionales suelen estar incompletos o ser observados de forma irregular.
Con respecto a la reconstrucción de las curvas, en la Tabla 5 es posible ver que el proceso de estimación tiene un mejor comportamiento cuando no se tiene error de medición, pues los valores error cuadrático medio relativo son bajos. No obstante, cuando se agrega error de medición, el método presenta mayor dificultad. Es posible evidenciar esto de una mejor forma en el escenario de los valores propios con decaimiento exponencial.
\[ \psi_m^{(2)}(t) = \sqrt{1-\alpha}\, f_m^{(2)}(t), \qquad t\in \mathcal{T}_2 := [-1,1]. \]
Para esta tercera simulación, se comparan el MFPCA calculado con el algoritmo FCP–TPA y otro método basado en B–splines en el cual se calcula el MFPCA con \(M_1=20\) imágenes propias y \(M_2=15\) funciones propias univariadas. En el caso de la estimación general de las imágenes y funciones propias, se realizó mediante productos tensoriales de \(K_1=10\times12\) funciones de bases de B–splines y \(K_2=15\) funciones, respectivamente.
Expansions <- function(SimData, SimData2, TrueFuns, M, uniExpansions){
TrueValues <- SimData$trueVals[1:M]
TrueFunctions <- TrueFuns
SimData3 <- SimData2
MFPCAObj <- MFPCA(SimData3,
M = M,
uniExpansions = uniExpansions)
FittedValues <- MFPCAObj$values
FittedFunctions <- MFPCAObj$functions
# Eigenvalue relative error
ErrValues <- (TrueValues - FittedValues) ^ 2 / (TrueValues ^ 2)
# Eigenfunction relative error
ErrFunsA <- norm(TrueFunctions - FittedFunctions, squared = TRUE)
ErrFunsB <- norm(TrueFunctions + FittedFunctions, squared = TRUE)
ErrFunsC <- norm(TrueFunctions, squared = TRUE)
ErrFuns <- pmin(ErrFunsA, ErrFunsB)
return(list(ErrValues = ErrValues,
ErrFuns = ErrFuns))
}Sim_Setting3 <- function(eValType = "exponential",
s2 = 0,
seed = NULL,
N = 250,
M,
Ms,
argvals = argvalsSetting3,
B){
uniExpansionsBasis <- list(list(type = "splines2D",
k = c(10, 12)),
list(type = "splines1D",
k = 15))
if(s2 > 0){
uniExpansionsBasis <- list(list(type = "splines2D",
k = c(10, 12)),
list(type = "splines1Dpen",
k = 15))
}
alphaRange = list(v = c(10^-5, 10^5),
w = c(10^-5, 10^5))
uniExpansionsTensor <- list(list(type = "FCP_TPA",
npc = 20,
alphaRange = alphaRange),
list(type = "uFPCA",
npc = 15))
ErrValues_Splines <- matrix(NA, nrow = B, ncol = M)
ErrFuns_Splines <- matrix(NA, nrow = B, ncol = M)
ErrValues_Tensor <- matrix(NA, nrow = B, ncol = M)
ErrFuns_Tensor <- matrix(NA, nrow = B, ncol = M)
set.seed(seed)
## MFPCA
pb <- txtProgressBar(min = 0, max = B, style = 3)
on.exit(close(pb), add = TRUE)
for (i in 1:B) {
SimData <- simMultiFunData(type = "weighted",
argvals = argvals,
M = Ms,
eFunType = list(c("Fourier", "Fourier"),
c("Poly")),
eValType = eValType,
N = N)
TrueFuns <- multiFunData(
funData(argvals = SimData$trueFuns[[1]]@argvals,
X = SimData$trueFuns[[1]]@X[1:M, , , drop = FALSE]),
funData(argvals = SimData$trueFuns[[2]]@argvals,
X = SimData$trueFuns[[2]]@X[1:M, , drop = FALSE])
)
if (s2 > 0) {
SimData2 <- addError(SimData$simData, sd = sqrt(s2))
} else {
SimData2 <- SimData$simData
}
Splines <- Expansions(SimData = SimData,
SimData2 = SimData2,
TrueFuns = TrueFuns,
M = M,
uniExpansions = uniExpansionsBasis)
Tensor <- Expansions(SimData = SimData,
SimData2 = SimData2,
TrueFuns = TrueFuns,
M = M,
uniExpansions = uniExpansionsTensor)
ErrValues_Splines[i, ] <- Splines$ErrValues
ErrFuns_Splines[i, ] <- Splines$ErrFuns
ErrValues_Tensor[i, ] <- Tensor$ErrValues
ErrFuns_Tensor[i, ] <- Tensor$ErrFuns
setTxtProgressBar(pb, i)
}
ErrValues = list(Splines = ErrValues_Splines,
Tensor = ErrValues_Tensor)
ErrFuns <- list(Splines = ErrFuns_Splines,
Tensor = ErrFuns_Tensor)
return(list(ErrValues = ErrValues,
ErrFuns = ErrFuns))
}A continuación se muestran una observación provenientes del dataset simulado, junto con las funciones propias que generan dichas observaciones mediante la descompsición de Karhúnen-Loève y los valores propios con decaimiento exponencial.
T1_1 <- seq(0, 1, length = 100)
T1_2 <- seq(0, 0.5, length = 50)
T2 <- seq(-1, 1, length = 200)
argvalsSetting3 <- list(list(T1_1, T1_2),
list(T2))
Ms <- list(c(5, 5), 25)
M <- 12
N <- 250
s2 <- 0.25
uniExpansionsBasis <- list(list(type = "splines2D",
k = c(10, 12)),
list(type = "splines1D",
k = 15))
alphaRange = list(v = c(10^-5, 10^5),
w = c(10^-5, 10^5))
uniExpansionsTensor <- list(list(type = "FCP_TPA",
npc = 20,
alphaRange = alphaRange),
list(type = "uFPCA",
npc = 15))
SimDataS3 <- simMultiFunData(type = "weighted",
argvals = argvalsSetting3,
M = Ms,
eFunType = list(c("Fourier", "Fourier"),
c("Poly")),
eValType = "exponential",
N = N)
SimDataS3s2 <- addError(SimDataS3$simData, sd = sqrt(s2))
Ahora, se muestra el procedimiento con el cual se realiza la simulación y se obtienen los resultados.
T1_1 <- seq(0, 1, length = 100)
T1_2 <- seq(0, 0.5, length = 50)
T2 <- seq(-1, 1, length = 200)
argvalsSetting3 <- list(list(T1_1, T1_2),
list(T2))
Ms <- list(c(5, 5), 25)
M <- 12
N <- 250
s2 <- 0.25
S3_exp_s0 <- Sim_Setting3(eValType = "exponential",
s2 = 0,
seed = 3,
N = 250,
M,
Ms,
argvals = argvalsSetting3,
B = 100)
S3_exp_s <- Sim_Setting3(eValType = "exponential",
s2 = s2,
seed = 3,
N = 250,
M,
Ms,
argvals = argvalsSetting3,
B = 100)En la Figura 49 se pueden ver los resultados obtenidos para esta tercer configuración.
Para la comparación de ambos métodos, se toman las primeras 12 componentes principales. Como se puede ver, en ambos escenarios (un decaimiento exponencial y lineal de los valores propios), el método de B–splines estima mucho mejor los valores e imágenes/funciones propias. En la presencia de error, los dos tienen resultados similares, pero el método de expansiones univariadas de B–splines se comporta mejor. Además, es posible ver que los últimos valores propios e imágenes/funciones propias son más difíciles de estimar con ambos enfoques.
En la quinta sección del paper se presenta una aplicación a un dataset que contiene curvas e imágenes relacionadas a un estudio ADNI (Alzheimer’s Disease Neuroimaging Initiative) que busca identificar biomarcadores que ayuden a la precisión de los diagnósticos de Alzheimer en una etapa temprana. Cada persona perteneciente al estudio tiene dos tipos de información:
Luego, por cada individuo se tiene asociado un vector \(\mathbf{X}_i=(X_i^{(1)}, X_i^{(2)})^t\) en el que la primera componente representa una trayetoria ADAS-Cog y la segunda una imágen cerebral FDG-PET.
Al intentar replicar el ejercicio realizado en el paper, se tuvo acceso a los datos de \(N=236\) participantes que tenían una imágen cerebral FDG-PET en la visita basal (una imagen cerebral al inicio del estudio) y almenos tres mediciones de ADAS-Cog durante el seguimiento. Debido a que las trayectorias no necesariamente estaban completas para todos los tiempos y todos los individuos, se tienen datos irregulares o sparse. Esto se muestra en la Figura 50.
Por otro lado, también se muestra un ejemplo de una de las imágenes cerebrales FDG-PET.
Para realizar el análisis de estos datos, primero se crean los objetos funcionales con los cuales se trabajará.
ADAS_fd <- funData(argvals = argvals_adas,
X = X_ADAS)
PET_fd <- funData(argvals = list((x_ind - 1) * 1.5,
(y_ind - 1) * 1.5),
X = PET_array_final)
names(ADAS_fd) <- rid_names
names(PET_fd) <- rid_names
ADNI_fd <- multiFunData(ADAS_fd, PET_fd)Además, es importante ser conscientes de que los tipos de datos son bastante diferentes con respecto a la escala y dominio. Por lo tanto, es importante es tomar ponderaciones o una de las dos partes podría dominar artificialmente el análisis. Para ello, se usa el producto escalar ponderado \[ \langle\langle f,g\rangle\rangle_w = \sum_{j=1}^{p} w_j \int_{T_j} f^{(j)}(t_j)g^{(j)}(t_j)\,dt_j \] para dos funciones \(f,g \in \mathcal{H}\). Los autores realizan la ponderación de la siguiente forma: \[ w_j = \left( \int_{T_j} \widehat{C}_{jj}(t_j,t_j)\,dt_j \right)^{-1} \] con \(j=1,2\), que correspondería a la variabilidad de cada uno de los elementos (trayectorias e imágenes cerebrales). Por lo tanto, los elementos reescalados \(\widetilde{X}^{(j)} = w_j^{1/2}X^{(j)}\) tienen varianza 1. A continuación se muestran los pesos:
ADAS_var_fd <- funData(
argvals = ADAS_fd@argvals,
X = array(apply(ADAS_fd@X, 2, var, na.rm = TRUE),
dim = c(1, length(argvals_adas)))
)
PET_var_fd <- funData(
argvals = PET_fd@argvals,
X = array(apply(PET_fd@X, 2:3, var, na.rm = TRUE),
dim = c(1, length(x_ind), length(y_ind)))
)
normScores <- integrate(ADAS_var_fd)
normPET <- integrate(PET_var_fd)
weights_ADNI <- c(1 / normScores, 1 / normPET)
weights_ADNI[1] 0.0001047925 0.0007438778De esta forma, se aplica la metodología propuesta en el paper. En el caso de las trayectorias, se realiza FPCA con \(M_1=3\) componentes principales que explican el 99% de la varianza. Para las imágenes, se usan productos tensoriales de \(20 \times 15\) B–splines.
ADNI_MFPCA <- MFPCA(
ADNI_fd,
M = 10,
uniExpansions = list(list(type = "uFPCA",
nbasis = 5,
npc = 3),
list(type = "splines2D",
k = c(20, 15))),
weights = weights_ADNI,
fit = FALSE,
approx.eigen = FALSE,
bootstrap = TRUE,
nBootstrap = 100,
bootstrapAlpha = c(0.05, 0.1),
verbose = TRUE
)| Componente | Valor Propio | Varianza explicada (%) | Varianza acumulada (%) |
|---|---|---|---|
| 1 | 0.6432 | 52.22 | 52.22 |
| 2 | 0.3815 | 30.98 | 83.20 |
| 3 | 0.0483 | 3.93 | 87.12 |
| 4 | 0.0378 | 3.07 | 90.19 |
| 5 | 0.0345 | 2.80 | 92.99 |
| 6 | 0.0257 | 2.08 | 95.07 |
| 7 | 0.0192 | 1.56 | 96.63 |
| 8 | 0.0167 | 1.35 | 97.99 |
| 9 | 0.0133 | 1.08 | 99.07 |
| 10 | 0.0114 | 0.93 | 100.00 |
Como se puede ver las dos primeras componentes ya explican más del 80% de la variabilidad ponderada, mientras que las cuatro primeras explican un poco más del 90%. Si nos concentramos en la primera componente que corresponde a un cambio en la media, entonces esta representa personas cuya función cognitiva empeora durante el seguimiento. Además, las personas asociadas positivamente con este componente tienden a mostrar menor actividad cerebral en regiones relacionadas con Alzheimer ya desde la visita basal.
Las imágenes de la primera función propia estimada \(\hat{\psi}_1\) se puede ver en la Figura 52. Es posible apreciar que la primera función propia presenta valores negativos en la componente ADAS-Cog durante todo el seguimiento. Esto quiere decir que sujetos con scores positivos se asocian con un menor deterioro cognitivo relativo (trayectorias ADAS-Cog por debajo de la media), mientras que sujetos con scores negativos se asocian con mayor deterioro cognitivo relativo (trayectorias por encima de la media). En cuanto a la componente FDG-PET, las regiones significativas positivas implican que los sujetos con scores negativos presentan menor señal FDG-PET relativa en dichas regiones.
Por otro lado, en la Figura 53 es posible ver que los participantes con diagnóstico AD (personas con Alzheimer diagnosticado) tienen scores más negativos, es decir, un deterioro cognitivo mayor. Por otro lado, las personas sanas (Normal) tienden a tener scores positivos, es decir, un menor deterioro cognitivo. Por último, aquellos participantes con deterioro cognitivo leve (MCI) se ubica en una posición intermedia.
En el caso de la estimación de la segunda función propia \(\hat{\psi}_2\) presente en la Figura 54, hay que tener en cuenta que las bandas de confianza de \(\hat{\psi}_1^{(1)}\) contienen al cero al 95% (y también 90%) de confianza. Es decir que no hay suficiente evidencia estadística para concluir que esta componente esté asociada con cambios importantes en ADAS-Cog (deterioro cognitivo). Por otro lado, en la segunda componente \(\hat{\psi}_1^{(2)}\) se puede ver que gran parte de la imágen propia tiene regiones significativas bastante extensas, sugiriendo que la segunda función propia está dominada principalmente por variación en las imágenes cerebrales FDG-PET.
Además, se presentan los scores estimados por género, asociados a la segunda componente principal en la Figura 55. De estos resultados es interesante ver la separación entre ambos grupos (Mujeres y Hombres) pues la mayoría de mujeres suele tener un score negativo, mientras que los hombres suelen tener scores positivos. En el paper, esto se asocia a las diferencias de registro o tamaño cerebral de los participantes en el estudio que se encuentra correlacionada con el género.
