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)load("Lambdas.RData")
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)
}
Uno de los supuestos para usar FPCA es que la media del proceso permanezca constante. Para verificar este supuesto, se asume que las observaciones funcionales \(\chi_i \in L^2\) son independientes. El objetivo es contrastar si la media de dichas funciones permanece constante en el tiempo \(i\). Por tanto, se define la hipótesis nula como:
\[ \begin{equation} H_0 : E [\chi_1] = E [\chi_2] = \dots = E [\chi_N] \end{equation} \] Y cada observación funcional puede representarse mediante el modelo:
\[ \begin{equation} \chi_i(t) = \mu(t) + Y_i(t), \quad E[Y_i(t)] = 0 \end{equation} \] Donde \(\mu(t)\) es la media poblacional.
Estadístico de Prueba
Para explicar la idea del procedimiento de prueba, se define, de acuerdo con (Horváth y Kokoszka 2012) y (Berkes et al. 2009)
\[ \hat{\mu}_k(t) = \frac{1}{k} \sum_{1 \leq i \leq k} X_i(t), \qquad \check{\mu}_k(t) = \frac{1}{N-k} \sum_{k < i \leq N} X_i(t). \]
Si la media es constante, la diferencia
\[ \Delta_k(t) = \hat{\mu}_k(t) - \check{\mu}_k(t) \]
es pequeña para todo \(1 \leq k < N\) y para todo \(t \in [0,1]\).
Sin embargo, \(\Delta_k(t)\) puede volverse grande debido a la variabilidad si \(k\) está cerca de \(1\) o de \(N\). Por esta razón, es usual trabajar con
\[ P_k(t) = \sum_{1 \leq i \leq k} X_i(t) - \frac{k}{N} \sum_{1 \leq i \leq N} X_i(t) = \frac{k(N-k)}{N} \left[ \hat{\mu}_k(t) - \check{\mu}_k(t) \right]. \]
Ya que el valor de \(P_k(t)\) no cambia si las funciones \(X_i(t)\) son reemplazadas por \(X_i(t)-\bar{X}_N(t)\). Al definir \(k=[Nx]\), con \(x \in (0,1)\), se obtiene
\[ \begin{equation} \int \left\{ \sum_{1 \leq i \leq Nx} X_i(t) - \frac{[Nx]}{N} \sum_{1 \leq i \leq N} X_i(t) \right\} \hat{v}_{\ell}(t)\,dt = \sum_{1 \leq i \leq Nx} \hat{\xi}_{\ell,i} - \frac{[Nx]}{N} \sum_{1 \leq i \leq N} \hat{\xi}_{\ell,i}. \end{equation} \] Lo que muestra que los scores pueden utilizarse para probar la constancia de la función de media.
El siguiente teorema puede emplearse para derivar varios estadísticos de prueba.
\[ \begin{equation} T_N(x) = \frac{1}{N} \sum_{\ell=1}^{d} \hat{\lambda}_{\ell}^{-1} \left( \sum_{1 \leq i \leq Nx} \hat{\xi}_{\ell,i} - x \sum_{1 \leq i \leq N} \hat{\xi}_{\ell,i} \right)^2 \end{equation} \] y sea \(B_1(\cdot), \ldots, B_d(\cdot)\) puentes brownianos estándar independientes.
Bajo \(H_0\),
\[ T_N(x) \overset{d}{\longrightarrow} \sum_{1 \leq \ell \leq d} B_{\ell}^{2}(x), \qquad (0 \leq x \leq 1), \]
El estadístico de prueba tipo Cramér-von Mises, denotado como \(S_{N,d}\), se define para \(d\) componentes principales como:
\[ \begin{equation} S_{N,d} = \frac{1}{N^2} \sum_{\ell=1}^{d} \frac{1}{\hat{\lambda}_\ell} \sum_{k=1}^{N} \left( \sum_{i=1}^k \hat{\xi}_{\ell,i} - \frac{k}{N} \sum_{i=1}^N \hat{\xi}_{\ell,i} \right)^2 \end{equation} \] Bajo la hipótesis nula, el estadístico \(S_{N,d}\) converge en distribución a la suma de cuadrados de puentes brownianos independientes:
\[ \begin{equation} S_{N,d} \xrightarrow{d} \int_0^1 \sum_{\ell=1}^d B_\ell^2(x) dx \end{equation} \]
La decisión estadística se toma comparando el valor calculado de \(S_{N,d}\) con los valores críticos. Se rechaza \(H_0\) con un nivel de significancia \(\alpha\) si: \[ \begin{equation} S_{N,d} > C_d(\alpha) \end{equation} \] Teniendo en cuenta que \(x \in [0, 1]\), el valor \(x^*\) que maximiza \(T_N(x)\), multiplicado por \(N\) corresponde a un estimador consistente del punto de cambio.
Este procedimiento se llevara a cabo a través la función fchange() del paquete fChange, donde es posible aplicar la metodología propuesta.
X.dfts <- dfts(fits[ , , 2])
N <- ncol(X.dfts$data)
ChangeDetectionSeg <- fchange(X.dfts,
method = "mean",
statistic = "Tn",
type = "segmentation",
critical = "simulation")| Punto de Cambio | p-valor |
|---|---|
| 134 | 0.000 |
| 184 | 0.004 |
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.
El Grupo_1 tiene 134 curvas.
El Grupo_2 tiene 50 curvas.
El Grupo_3 tiene 84 curvas.
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 3.
X <- t(fits[, , 2])
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 el proceso, se encuentran las estimaciones de las funciones propias asociadas al operador de covarianza del proceso funcional.
Las tres primeras funciones propias se pueden visualizar en la Figura 4.
[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 5, Figura 6 y Figura 7.
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.
En las gráficas anteriores es posible ver cuánto porcentaje de variabilidad es explicado por las tres funciones propias estimadas. Esto también se muestra a través de un scree plot en la Figura 8.
De esta forma, es posible ver que 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 el 96.33 de la variabilidad, aproximadamente. Por lo tanto, para cada observación funcional, se tiene su aproximación \[ X_i(t) \approx \hat \mu(t) + \sum_{k=1}^K \hat {\xi}_{ik} \hat \nu _ k(t) \] con \(i = 1,\ldots,268\) y \(K=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, "+") 
A continuación, en la Tabla 2 se muestra el error cuadrático medio de predicción para cada observación.
| MPSE | |
|---|---|
| Primera Observación | 1.2944 |
| Segunda Observación | 3.8829 |
| Tercera Observación | 7.1369 |
A continuación se presentan las realizaciones del primer score en la Figura 12.
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 Figura 14 se presenta la estimación de la función de covarianza y las funciones propias.
cov.fd <- var.fd(X.fd$fd)
cov.grid <- seq(min(t), max(t), length = 20)
cov.mat <- eval.bifd(cov.grid, cov.grid, cov.fd)
Como es posible visualizar, 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 Figura 15 se puede ver la función de covarianza junto con su aproximación y en la Figura 16 se ven sus residuales.
eigenfn2 <- eval.fd(cov.grid, pcaObj$harmonics)
NU <- as.matrix(eigenfn2)
D <- diag(eigenval)
mercer <- NU %*% D %*% t(NU)
cov.diff <- mercer - cov.mat
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.
Ya que se identificaron puntos de cambio, la construcción de las cartas de control se llevará a cabo mediante dos estrategias de partición del conjunto de datos, ambas basadas en la selección de observaciones que comparten una misma 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 <- sample(1:n_total, size = floor(r * n_total), replace = FALSE)
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
))
}En la primera opción, se utilizará como Fase I el subconjunto con mayor número de observaciones bajo un comportamiento estable, el cual en este caso corresponde a las primeras 134 observaciones. A partir de este conjunto se estimarán los límites de control mediante bootstrap. Posteriormente, el conjunto restante de observaciones será utilizado como Fase II, con el fin de evaluar el comportamiento del proceso bajo los límites previamente construidos.
resT2Q_1=clean_phaseI_T2Q (data=wvl, p=2, argvals=t, Lambdas=Lambdas, Basis=Basis,B = 500,r=1,cp=1:134)De las 134 observaciones inicialmente consideradas en la Fase I, 59 permanecieron bajo control y fueron utilizadas para la estimación final de los límites de control. A partir de estas observaciones, se obtuvieron como límites 5.682479 para el estadístico \(T^2\) y 341.1768 para el estadístico \(Q\).
| Observación | Valor T² |
|---|---|
| 151 | 6.112 |
| 166 | 7.328 |
| 178 | 7.075 |
| 190 | 22.086 |
| 199 | 6.659 |
| 231 | 5.904 |
| 268 | 26.473 |
El porcentaje de observaciones fuera de control para el estadístico Q es de 47.76%.En la segunda opción, se considerará una partición alternativa del conjunto de datos, tomando como Fase I las últimas observaciones que también presentan 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. Al igual que en el caso anterior, los límites de control serán obtenidos 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, 50 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.571344 para el estadístico \(T^2\) y 300.8219 para el estadístico \(Q\).
| Observación real | Valor T² | Valor Q | Fuera por T² | Fuera por Q |
|---|---|---|---|---|
| 185 | 1.571 | 1452.490 | No | Sí |
| 187 | 2.731 | 608.605 | No | Sí |
| 212 | 1.327 | 343.029 | No | Sí |
| 221 | 0.346 | 418.832 | No | Sí |
| 228 | 1.069 | 1050.650 | No | Sí |
| 230 | 0.511 | 736.786 | No | Sí |
| 231 | 5.904 | 223.268 | Sí | No |
| 266 | 0.644 | 475.134 | No | 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 1) y una fase de monitoreo (fase 2).
En la fase 1 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 1, mientras que el 20% restante se destinará a la fase 2.
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 <- sample(
1:n_total,
size = floor(r * n_total),
replace = FALSE
)
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, 39 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.103766
| Observación real | Valor DFM |
|---|---|
| 231 | 2.43 |
El supuesto esencial de esta metodología es la existencia en cada tiempo \(t\) de unas direcciones principales que pueden cambiar a lo largo del tiempo pero de manera suave. Por ejemplo, la temperatura del agua a una profundidad fija y la temperatura ambiente son dos variables que están asociadas, pero es posible que el nivel de asociación cambie a lo largo del día; el supuesto está en que si bien se contemplan cambios en los niveles de asociación, esos cambios deben ser suaves, explícitamente deben ser modelados por funciones diferenciables como se expuso en la construcción del espacio de trabajo.
El producto interno en traspone las mediciones de distancia a un espacio en el que no hay asociación.
-Existe un operador \(\mathcal{F}\) (Operador de covarianza) invertible y un proceso funcional \(p-\)variado de media 0 \[ \begin{equation*} \widetilde{\mathbf{Y}}_i = \left( \widetilde{Y}_{i,1}, \widetilde{Y}_{i,2}, \dots, \widetilde{Y}_{i,p} \right), \end{equation*} \] Las observaciones son realizaciones de un proceso funcional \(\mathbf{X}_i\) con \[ \begin{equation} \mathbf{X}_i = \boldsymbol{\mu}_{\gamma} + \mathbf{Y}_i \end{equation} \]
\(\mathbf{Y}_i\) es un proceso de media 0, \(E(\mathbf{X}_i) = \boldsymbol{\mu}_{\gamma}\) con \(\gamma = 1, 2\), y se tiene
\[ \begin{equation} \mathbf{Y}_i = \mathbf{X}_i - \boldsymbol{\mu}_{\gamma} \end{equation} \] donde
\[
\begin{equation}
\mathbf{Y}_i = \mathcal{F} \widetilde{\mathbf{Y}}_i, \quad \mathbf{X}_i = \mathcal{F} \widetilde{\mathbf{X}}_i, \quad \boldsymbol{\mu}_{\gamma} = \mathcal{F} \widetilde{\boldsymbol{\mu}}_{\gamma}
\end{equation}
\]
Si existe un punto de cambio en \(c\) igual a parte entera de \(N\theta\) con \(\theta \in [0,1]\) las observaciones provienen del proceso \(\mathbf{X}_i\) en el que desde 1 hasta el valor entero de \(c\) se tiene \(\boldsymbol{\mu}_{\gamma} = \boldsymbol{\mu}_1\) y después \(\boldsymbol{\mu}_{\gamma} = \boldsymbol{\mu}_2\). Si no hay punto de cambio \(\boldsymbol{\mu}_1 = \boldsymbol{\mu}_2\).
Para cada \(k\) las realizaciones \(Y_{i,k}\) provienen del proceso \(Y_k\) que es una componente del proceso funcional \(p-\)variado \(\mathbf{Y}_i\). Siendo este último la transformación de un proceso \(\widetilde{\mathbf{Y}}_i\) como se indicó en el anterior ítem.
Considérese la siguiente muestra de tamaño \(N\):
\[ \begin{equation} \mathbf{X}_1, \mathbf{X}_2, \dots, \mathbf{X}_N \end{equation} \] con \[ \begin{equation} \mathbf{X}_i = (X_{i,1}, X_{i,2}, \dots, X_{i,p}), \end{equation} \] para \(i = 1, 2, \dots, N\).
Las siguientes definiciones son análogas a las expuestas en
\[ \begin{equation} \hat{\mu}_{c,k}(t) = \frac{1}{c} \sum_{i=1}^{c} X_{i,k}(t), \quad \hat{\boldsymbol{\mu}}_c = (\hat{\mu}_{c,1}(t), \hat{\mu}_{c,2}(t), \dots, \hat{\mu}_{c,p}(t)) \tag{4-51} \end{equation} \] \[ \begin{equation} \check{\mu}_{c,k}(t) = \frac{1}{N-c} \sum_{i=c+1}^{N} X_{i,k}(t), \quad \check{\boldsymbol{\mu}}_c(t) = (\check{\mu}_{c,1}(t), \check{\mu}_{c,2}(t), \dots, \check{\mu}_{c,p}(t)) \tag{4-52} \end{equation} \] Ahora se considera la diferencia entre los estimadores de las medias \(\hat{\boldsymbol{\mu}}_c\) y \(\check{\boldsymbol{\mu}}_c\) multiplicado por el factor de tipo parabólico.
\[ \begin{equation} P_c(t) = \frac{c(N-c)}{N} (\hat{\boldsymbol{\mu}}_c(t) - \check{\boldsymbol{\mu}}_c(t)) \tag{4-53} \end{equation} \]
el reemplazo de \(\hat{\boldsymbol{\mu}}_c\) y \(\check{\boldsymbol{\mu}}_c\) en la anterior expresión, conlleva a \[ \begin{equation} P_c(t) = \sum_{1 \le i \le Nx} (\mathbf{X}_i(t) - \overline{\mathbf{X}}_N(t)) - \frac{c}{N} \sum_{i=1}^{N} (\mathbf{X}_i(t) - \overline{\mathbf{X}}_N(t)) \tag{4-54} \end{equation} \] donde \[ \begin{equation} \overline{\mathbf{X}}_N = \frac{1}{N} \sum_{i=1}^{N} \mathbf{X}_i \end{equation} \]
Se define \(P_t(x)\) como función de \(x \in (0,1)\) mediante \[ \begin{equation} P_t(x) = \sum_{1 \le i \le Nx} (\mathbf{X}_i(t) - \overline{\mathbf{X}}_N(t)) - x \sum_{i=1}^{N} (\mathbf{X}_i(t) - \overline{\mathbf{X}}_N(t)) \end{equation} \]
Para \(\mathbf{X}_i\), se tiene que \(P_t(x) \in \mathcal{F}\mathbf{H}_D\).
Se define
\[ \begin{equation*} T_N(x) := \frac{1}{N} \| P_t(x) \|^2_{\hat{K}} \end{equation*} \] Sea
\[ \begin{equation} \hat{\eta}_{i,m,k} = \langle \widetilde{\mathbf{X}}_i - \overline{\widetilde{\mathbf{X}}}_N, \hat{\mathbf{e}}_{m,k} \rangle \tag{4-55} \end{equation} \]
Finalmente, la expresión para la norma al cuadrado se puede escribir como:
\[ \begin{equation} \| P_t(x) \|^2_{\hat{K}} = \sum_{k=1}^{p} \sum_{m=1}^{d_k} \left[ \frac{1}{\hat{\lambda}_{m,k}} \left( \sum_{1 \le i \le Nx} \hat{\eta}_{i,m,k} - x \sum_{i=1}^{N} \hat{\eta}_{i,m,k} \right)^2 \right] \tag{4-59} \end{equation} \]
Para determinar las componentes principales funcionales multivariadas, se seguirá el enfoque propuesto por Clara Happ, 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)ChangeMultiData <- list(procesos = aperm(fits, c(1, 3, 2)),
nombres = paste0("var", 1:dim(fits)[3]),
tiempos = grid)
ChangeDetectionMulti <- detectar(datos = ChangeMultiData,
no_tiempos = "dim(ChangeMultiData$procesos)[1]",
no_bases = "40",
com_prin = "rep(3, 7)")
load("MFPCAObj.RData")
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, C., & Greven, S. 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.
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.
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 <- 8SimDataS1 <- 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 26.
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, 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().
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))
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.
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 30.
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))
}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)
PETMETA <- read.csv("PETMETA_ADNI1.csv")
ADNI_Meta <- read.csv("ADNI_Meta.csv")
PET_IDs <- unique(PETMETA$RID)
meses <- data.frame(
VISCODE = c("bl","m06","m12","m18","m24","m36","m48","m60","m72","m84","m96"),
Month = c(0, 6, 12, 18, 24, 36, 48, 60, 72, 84, 96)
)
ADNI_plot <- ADNI_Meta %>%
filter(RID %in% PET_IDs) %>%
filter(!is.na(ADASCog)) %>%
filter(VISCODE %in% meses$VISCODE) %>%
left_join(meses, by = "VISCODE") %>%
group_by(RID) %>%
filter(n() >= 3) %>%
ungroup()
# Conteo de observaciones por visita
conteos <- ADNI_plot %>%
group_by(Month) %>%
summarise(n = n_distinct(RID), .groups = "drop")
ggplot(ADNI_plot, aes(x = Month, y = ADASCog, group = RID)) +
geom_line(color = "grey75", linewidth = 0.3) +
geom_point(color = "grey65", size = 0.8) +
geom_text(
data = conteos,
aes(x = Month, y = min(ADNI_plot$ADASCog, na.rm = TRUE) - 2, label = n),
inherit.aes = FALSE,
size = 3
) +
scale_x_continuous(breaks = meses$Month) +
labs(
x = "Months since baseline",
y = "ADAS-Cog"
) +
theme_classic()
load("Application_Data.RData")
rid_names <- as.character(rownames(X_ADAS))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_namesADNI_fd <- multiFunData(ADAS_fd, PET_fd)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.0007438778ADNI_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
)load("ADNI_MFPCA.RData")
plot(ADNI_MFPCA$functions, obs = 1)
plot(ADNI_MFPCA$functions, obs = 2)
scores_df <- as.data.frame(ADNI_MFPCA$scores)
colnames(scores_df) <- paste0("Score", seq_len(ncol(scores_df)))
scores_df$RID <- as.integer(rownames(X_ADAS))
head(scores_df) Score1 Score2 Score3 Score4 Score5 Score6
2 0.2060327 1.2704571 0.40815676 0.2875352 0.01131930 0.07970470
3 -0.7460181 1.2005240 0.18038238 0.2076549 0.00152945 -0.12092085
5 0.4422881 0.1102053 0.14320636 -0.1439783 -0.21832905 -0.05835144
8 0.3589889 -0.9122447 0.08578808 0.1698420 0.16218625 0.09495979
10 -0.5518996 -0.3365288 0.03999698 -0.1749023 0.03146703 -0.11270949
16 0.3683285 -0.7030324 -0.06157218 -0.1362819 -0.23220698 -0.05307746
Score7 Score8 Score9 Score10 RID
2 0.36805824 -0.26262578 0.05379191 -0.23931872 2
3 0.06106803 0.04713955 -0.09307684 0.05219218 3
5 0.16581187 0.06351632 -0.04218105 -0.04882512 5
8 -0.14352340 -0.03764256 -0.20039436 -0.07834631 8
10 0.03647830 0.15682456 -0.15248548 -0.04518432 10
16 0.19615744 0.03448786 0.09791013 0.12799415 16meta_subject <- ADNI_Meta %>%
filter(VISCODE == "bl") %>%
distinct(RID, .keep_all = TRUE) %>%
select(RID, DIAGNOSIS, GENDER)
scores_meta <- scores_df %>%
left_join(meta_subject, by = "RID")
scores_long_diag <- scores_meta %>%
select(RID, DIAGNOSIS, starts_with("Score")) %>%
pivot_longer(
cols = starts_with("Score"),
names_to = "Component",
values_to = "Score"
)
scores_long_gender <- scores_meta %>%
select(RID, GENDER, starts_with("Score")) %>%
pivot_longer(
cols = starts_with("Score"),
names_to = "Component",
values_to = "Score"
)ggplot(
scores_long_diag %>% filter(Component %in%"Score1"),
aes(x = DIAGNOSIS, y = Score)
) +
geom_boxplot() +
facet_wrap(~ Component, scales = "free_y") +
theme_classic() +
labs(
x = "Diagnosis",
y = "MFPCA score",
title = "Boxplots of MFPCA scores by diagnosis"
)
ggplot(
scores_long_gender %>% filter(Component %in% "Score2"),
aes(x = GENDER, y = Score, group = GENDER)
) +
geom_boxplot() +
facet_wrap(~ Component, scales = "free_y") +
theme_classic() +
labs(
x = "Gender",
y = "MFPCA score",
title = "Boxplots of MFPCA scores by gender"
)