FDA REGRESION

Lectura del conjunto de datos

library(fda)
## Cargando paquete requerido: splines
## Cargando paquete requerido: fds
## Cargando paquete requerido: rainbow
## Cargando paquete requerido: MASS
## Cargando paquete requerido: pcaPP
## Cargando paquete requerido: RCurl
## Cargando paquete requerido: deSolve
## 
## Adjuntando el paquete: 'fda'
## The following object is masked from 'package:graphics':
## 
##     matplot
P<-read.csv("caso_5_split.csv")
head(P)

Se dispone de un conjunto de datos, denominado A, que engloba diversas mediciones (d40, d47, …, d82) efectuadas en intervalos de tiempo uniformes. Dichas mediciones están vinculadas con distintos tratamientos de fungicida y clones. Además, el conjunto de datos contempla dos categorías de producción: Comercial y No Comercial.

El objetivo es examinar la producción (tanto Comercial como No Comercial) en relación con el tiempo y determinar cómo diferentes tratamientos (en este caso, el fungicida 0) influyen en esta producción.

Protocolo

Procesamiento de Datos:

Se instauró una secuencia temporal que simboliza los intervalos regulares en los que se efectuaron las mediciones. Posteriormente, se aplicó un filtro al conjunto de datos A para extraer exclusivamente aquellos registros donde el fungicida es de tipo 0 (almacenados en F0). Seguidamente, se seleccionaron únicamente las columnas de mediciones de F0 para el análisis (almacenados en M).

Constitución de Bases de Fourier:

Se empleó la función create.fourier.basis para instaurar una base de Fourier, que es beneficiosa para modelar ciclos o patrones periódicos en los datos. La base de Fourier se instauró para el rango de tiempo de 40 a 82 con 5 bases (nbasis = 5).

Transformación de Datos a Funciones:

Con la función Data2fd, se transmutaron los datos en un objeto funcional fdRegre utilizando la base de Fourier que se instauró previamente.

Modelo de Regresión Funcional: Posteriormente, se instauró un modelo de regresión funcional fdModel utilizando la función fRegress, donde la producción total (la suma de Comercial y NoComercial) se modeló en función de las mediciones transmutadas en funciones (fdRegre).

Predicción: Se efectuaron predicciones con el modelo y se fusionaron con los datos de producción reales para comparar (B). Se exhiben los primeros 10.

tiempo<-seq(40,82,7)
F0<-subset(P,fungicida==0)
M<-F0[,6:12]
base <- create.fourier.basis(rangeval = c(40,82),nbasis = 5)
fdRegre<-Data2fd(argvals = tiempo, y=t(M),basisobj = base)
produccion <- as.numeric(F0$comercial + F0$NoComercial)
fdModel<-fRegress(produccion ~ fdRegre)
prediccion <- predict(fdModel)
cbind(produccion,prediccion)->B
plot(fdRegre, main = NULL, xlab = "Tiempo (días)", ylab = "Producción")

## [1] "done"
head(B,10)
##         produccion          
## reps 1       19.80 20.301410
## reps 2       30.00 30.513930
## reps 3       22.00 23.494978
## reps 4       40.10 35.601804
## reps 5       25.60 23.349681
## reps 6       20.40 24.138143
## reps 7        0.25  3.307956
## reps 8       19.70 18.041935
## reps 9       30.30 29.855791
## reps 10      27.60 29.995097

Base B-Spline

basebspline <- create.bspline.basis(rangeval = tiempo, norder = 9)
fdRegre2<-Data2fd(argvals = tiempo, y=t(M),basisobj = basebspline)
produccion <- as.numeric(F0$comercial + F0$NoComercial)
fdModel2<-fRegress(produccion ~ fdRegre2)
prediccion2 <- predict(fdModel2)
cbind(produccion,prediccion2)->B2
head(B2,10)
##         produccion          
## reps 1       19.80 21.973845
## reps 2       30.00 29.792913
## reps 3       22.00 22.747664
## reps 4       40.10 39.579246
## reps 5       25.60 24.390577
## reps 6       20.40 20.245815
## reps 7        0.25  1.773437
## reps 8       19.70 18.988075
## reps 9       30.30 30.242053
## reps 10      27.60 28.617930
dim(B)
## [1] 78  2
plot(1:78,produccion,col=4,type="b",axes=F,xlab="Tiempo (días)",ylab="Producción" ,
ylim=c(0, 65))
axis(1)
axis(2)
lines(1:78,prediccion,col=2)
lines(1:78,prediccion2,col=3)
legend("topright", legend=c("Producción", "Producción Fourier", "Producción BSpline"),
col=c(4, 2, 3), lty=1, pch=c(NA, NA, NA), bty="n")

Las mediciones auténticas de producción (ilustradas en azul) y las estimaciones del modelo (representadas en rojo) se extienden a lo largo de 78 unidades temporales. Los puntos azules, que simbolizan la producción real, exhiben una variabilidad notable, mientras que la línea roja, que denota las predicciones, sigue una tendencia más uniforme. La cercanía entre ambas series en diversos puntos insinúa que el modelo posee un grado de precisión considerable, aunque existen áreas donde las estimaciones no se alinean exactamente con los valores auténticos.

Metricas del modelo

O<-as.numeric(fdModel$yvec)
P<-as.numeric(fdModel$yhatfdobj)
n<-length(O)
p<-base$nbasis
e<-O-P
SSE <- sum(e^2)
RMSE <- sqrt(SSE/n)
SSY <- sum((O - mean(O))^2)
R2.fourier <- (SSY-SSE)/SSY
gl1 <- p-1 
gl2 <- n-p 
Fratio <- ((SSY-SSE)/gl1)/(SSE/gl2)
pvalue<-1-pf(Fratio,gl1,gl2)
R2.fourier
## [1] 0.9058508
plot(P,P,type="l")
points(P,O,pch=19,col="green")

Se observa una relació n lineal entre dos conjuntos de datos, una comparació n de valores predichos “P” contra sı́ mismos para formar una lı́nea de referencia. Los puntos verdes, que representan los valores observados “O”, están dispersos alrededor de esta lı́nea, indicando la precisió n de las predicciones. La cercanı́a de los puntos verdes a la lı́nea muestra que las predicciones son bastante precisas, aunque algunos puntos se desvı́an, lo que indica errores de prediccion.

BSpline

O2<-as.numeric(fdModel2$yvec)
P2<-as.numeric(fdModel2$yhatfdobj)
n2<-length(O2)
p2<-basebspline$nbasis
e2<-O2-P2
SSE <- sum(e2^2)
RMSE <- sqrt(SSE/n2)
SSY <- sum((O2 - mean(O2))^2)
R2.spline <- (SSY-SSE)/SSY
gl1 <- p2-1 
gl2 <- n2-p2
Fratio <- ((SSY-SSE)/gl1)/(SSE/gl2)
pvalue<-1-pf(Fratio,gl1,gl2)
R2.spline
## [1] 0.9820393

El coeficiente de determinación, conocido como R², ostenta un valor de 0.9820393, lo cual indica que aproximadamente el 98.20% de la variación observada en los datos de producción puede ser atribuida al modelo de regresión funcional empleando una base Spline. Esto implica que la capacidad explicativa del modelo respecto a la variabilidad de los datos es sumamente alta, reflejando una robusta correlación entre las variables independientes y la variable dependiente analizada.

plot(P2,P2,type="l")
points(P2,O2,pch=19,col="blue")

Comparacion del Coeficiente de Dterminacion

cotejo <- data.frame(
R2Fourier = R2.fourier ,
R2.spline = R2.spline
)
colnames(cotejo) <- c("R2 Fourier", "R2 Spline")
knitr::kable(cotejo)
R2 Fourier R2 Spline
0.9058508 0.9820393

Valores puntuales de la media funcional por cada mes

media_f<-mean.fd(fdRegre)
mp<-eval.fd(tiempo,media_f)
plot(media_f)
## [1] "done"
points(tiempo,mp,pch=19,col="orange")

La representación gráfica expone la media funcional (mf) de la serie temporal, delineada por una trayectoria curvilínea suavizada, y los valores individuales de dicha media (mp) correspondientes a cada mes, señalizados mediante puntos de tonalidad azulada. La curva se inaugura marginalmente por encima del valor 10, exhibe una declinación tenue entre los intervalos temporales de 45 y 50, asciende hasta un ápice que sobrepasa el valor 20 en torno al instante temporal 75 y posteriormente decrece con sutileza. Esta modalidad de visualización facilita la detección de tendencias o patrones subyacentes en los datos funcionales, tales como las fases de incremento o reducción en la media, proporcionando así una herramienta analítica valiosa para la elucidación del comportamiento global de los datos a través del espectro temporal.

Predicciones Nuevas

constList <- function(model,periodos){
basis <- model$xfdlist$const$basis
coefs <-matrix(rep(1,periodos),c(1,periodos))
dimnames(coefs) <- list("time",paste("reps",1:periodos))
reps <- paste("reps",1:periodos)
fdnames <- list(args="time",reps=reps,funs="values")
const <- list(coefs=coefs,basis=basis,fdnames=fdnames)
class(const) <- "fd"
return(const)
}
parcela_1 <-c(6, 4, 5, 2, 8, 6, 5)
parcela_2 <-c(5, 9, 7, 5, 4, 3, 4)
nuevas<-rbind(parcela_1,parcela_2)
nuevas.fd<-Data2fd(argvals=tiempo,y=t(nuevas),basisobj=base)
beta<-constList(fdModel,2)
yhat<-predict(fdModel,newdata = list(beta,nuevas.fd))
round(yhat,1)
##      [,1]
## [1,] 15.2
## [2,] 13.3
sum(parcela_1)
## [1] 36
sum(parcela_2)
## [1] 37
e_1 <- sum(parcela_1) - round(yhat,1)[1,]
e_2 <- sum(parcela_2) - round(yhat,1)[2,]
d_1 <- e_1/sum(parcela_1)
d_2 <- e_2/sum(parcela_2)
d_d <- (d_1 + d_2)
d_d
## [1] 1.218318