Análisis de datos funcionales (ADF) - Modelos de regresión

Curso de Estadística Espacial

Universidad Nacional de Colombia - Sede Bogotá

Modelos lineales funcionales (Regresión Funcional)

Para explicar la variabilidad de una determinada variable con respecto a otras explicativas, consideradas como covariables, el análisis de varianza y la regresión lineal son los procedimientos que generalmente se utilizan. En este sentido el modelo de regresión funcional es la extensión natural del modelo de regresión usual al caso en el cual se cuenta con una variable respuesta funcional o con covariables funcionales (Aristizabal, 2011).

Regresión lineal funcional con respuesta funcional (ANOVA Funcional)

En términos formales (Ramsay & Silverman, 2005), se asume que se cuenta con \(G\) “tratamientos” cada uno con un número \(N_g\) de sujetos. El modelo para la \(m\)-ésima función (curva) en el \(g\)-ésimo grupo (\(y_{mg}(t))\), esta dada por:\[ y_{mg}(t) = \mu(t) + \alpha_g(t) + \varepsilon_{mg}(t)\] dónde la función \(\mu\) es la media general, \(\alpha_g\) representa la función media para cada ``tratamiento’’ y \(\varepsilon_{mg}\) es la función de error en cada caso. La tarea entonces, es establecer la matriz diseño para lograr estimar los parámetros funcionales \(\mu\) y \(\alpha_g\), bajo la condición que \(\sum_g \alpha_g(t) = 0\) para todo \(t\) con el fin de garantizar la estimabilidad de los mismos.

En términos matriciales el modelo queda determinado como: \[y_{mg}(t) = \sum_{j=1}^{(G+1)} x_{(mg)j}\beta_j(t) + \varepsilon_{(mg)}(t) \]\[y(t) = X(t)\beta(t) + \varepsilon(t)\]

Cabe destacar que la matriz diseño \(X\) tiene la misma estructura que en el caso multivariado o univariado, la diferencia obedece a que el vector de parámetros \(\beta(t)\) y las predicciones \(X\beta(t)\), son vectores de funciones en vez de vectores de números.

Estimación - Mínimos cuadrados ordinarios

Para realizar el ajuste a través del criterio de mínimos cuadrados ordinarios, se debe escoger \(\beta(t)\) que minimize la suma de cuadrados residual, así, para extender el principio de mínimos cuadrados al caso funcional, se reinterpreta la suma de cuadrados residual (\(y_i(t) - X_i\beta(t)\)) como una nueva función, y entonces por el criterio de mínimos cuadrados ordinarios se debe minimizar:\[SSE(\beta) = \sum_g ^{G} \sum_m^{N_g} \int (y_{mg}(t) - \sum_j^q x_{(mg)j}\beta_j(t))^2dt = \sum_i^n \|y_i(t) - x_i^t\beta(t) \|^2\] Minimizando \(SSE(\beta)\) sujeto a la condición de \(\sum_{j = 2} ^{G+1}\beta_j(t) = 0\) se obtiene la estimación de los parámetros funcionales \(\hat{\beta}(t)\).

Al igual que en el modelo lineal multivariado, la fuente primaria de información para investigar la importancia de los ``tratamientos’’, es la función de suma de cuadrados:\[ SSE(t) = \sum_{mg} (y_{mg}(t) - X_{mg}\hat{\beta}(t))^2\] Está función es comparada con la función de suma de cuadrados de los errores obtenida al utilizar solo la media general \(\hat{\mu}\) en el modelo:\[SSY(t) = \sum_{mg} (y_{mg}(t) - \hat{\mu}(t))^2\] Así un camino para realizar la comparación, se establece utilizando la función del cuadrado de correlación múltiple:\[RSQ(t) = \frac{(SSY(t) - SSE(t))}{SSY(t)}\] o calcular las funciones análogas a las de una tabla ANOVA en el caso univariado. Por ejemplo la función de cuadrado medio del error:\[MSE = \frac{SSE}{df(error)} = \frac{SSE}{N-G}\] Análogamente, la función de cuadrado medio de la regresión es la diferencia entre la suma de cuadrados total (\(SSY\)) y la suma de cuadrados del error (\(SSE\)), dividida por la diferencia entre los grados de libertad del error para los dos modelos (grados de libertad de la regresión):\[MSR(t) = \frac{SSY(t) - SSE(t)}{df(regresion)}\] Finalmente se puede construir la función \(F\) como:\[ F = \frac{MSR}{MSE}\]

Básicamente, la mayoría de los estadísticos del análisis de varianza univariado son aplicables al problema funcional, esencialmente, porque el problema de análisis de varianza funcional, se puede considerar como un análisis de varianza univariado (ANOVA) para cada valor específico \(t\) del dominio de la función. Sin embargo, bajo este enfoque la prueba \(F\) pierde algunas de sus propiedades, debido a que realizar una prueba \(F\) en cada tiempo \(t\) con un nivel de significancia determinado no implica el mismo nivel de significancia en una prueba conjunta.

Pruebas de hipótesis - Prueba de Shen & Faraway

En el 2004 Shen y Faraway (Shen & Faraway, 2004), propusieron una prueba tipo \(\textbf{F}\) para modelos lineales con respuesta funcional. Esta, es una extensión de la prueba \(\textbf{F}\) multivariada al caso en que la cantidad de mediciones crece rápidamente y puede ser utilizada para la comparación de dos modelos anidados cualesquiera.

Sin pérdida de generalidad, se considera la comparación de dos modelos lineales \(\omega\) y \(\Omega\), dónde \(dim(\Omega) = p\) y \(dim(\omega) = q\), \(p>q\). El modelo \(\omega\) resulta de una restricción lineal de los parámetros de \(\Omega\). Así:\[ H_0: Y(t) = X_1\alpha_1(t) + \epsilon(t) \hspace{5pt}\text{versus} \hspace{5pt} H_a: Y(t) = X_1\alpha_1(t) + X_2\alpha_2(t) +\epsilon(t) \]

Sea \[\textbf{F} = \frac{(rss_\omega - rss_\Omega)/(p-q)}{rss_\Omega/(n-p)} \approx \frac{traza(\hat{\sum}^\omega - \hat{\sum}^\Omega)/(p-q)}{traza(\hat{\sum}^\Omega)/(n-p)}\] con \(rss = \sum_{i=1}^n \int_t (y_i(t) - \hat{y}_i(t))^2 dt.\)

Asumiendo que el proceso de ruido \(\epsilon_i(t)\) es gaussiano e independiente y con función de covarianza continua \(r(s,t)\) sobre un intervalo cerrado \(\tau\), entonces bajo el sistema de hipótesis descrito, la estadística se distribuye como \[\frac{\sum_{i=1}^\infty r_i \chi_{p-q}^2/(p-q)}{\sum_{i=1}^\infty r_i \chi_{n-p}^2/(n-p)}\] dónde \(r_i\) es el \(i\)-ésimo valor propio ordenado de la función de covarianza \(r(s,t)\) y las variables \(\chi^2\) son independientes.

En su artículo los autores demuestran que la estadística tiene una distribución denominada distribución funcional con coeficientes \({r_i, i =1,2,...,\infty}\) y grados de libertad \((p-q, n-p)\). Así mismo, realizan la aproximación de esta distribución a la distribución \(F\) de Fisher con grados de libertad \(f_1\) y \(f_2\):\[f_1 = \frac{(\sum_{i=1}^\infty r_i)^2}{\sum_{i=1}^\infty r_i^2}(p-q); f_2 = \frac{(\sum_{i=1}^\infty r_i)^2}{\sum_{i=1}^\infty r_i^2}(n-p)\] En la práctica, solo se observa \(y_i(t)\) sobre una grilla de puntos \(t_j\), \(j = 1,2,...,M\), así, el factor de ajuste \(\frac{(\sum_{i=1}^\infty r_i)^2}{\sum_{i=1}^\infty r_i^2}\) puede ser estimado por \(\frac{traza(E)^2}{traza(E^2)}\), dónde \(E = \hat{\sum}^\Omega\) es la matriz de covarianza empírica del ruido calculada con el modelo completo \(\Omega\).

Ejemplo

Experimento

La situación experimental llevada a cabo y que dió origen a la base de datos con la cual se realizan las aplicaciones, parte de que a un paciente se le mostraba un conjunto de palabras con carga emocional positiva (abrazo, paz, regalo), negativa (guerra, enfermedad, muerte), o neutra (número, roca), y también una serie de no palabras. Su tarea consistía en responder mediante un botón adyacente cuando aparecía únicamente una palabra. El número de pacientes total en el experimento fue de 48, de los cuales 30 son hombres y los restantes 18 mujeres.

La presentación de los estímulos y la realización de la tarea se hace en sincronía con el registro del electroencefalograma (EEG) con lo cual se obtiene la medición de la actividad del cerebro medida como la diferencia entre el voltaje máximo y la línea base preestimular (electrodo determinado) para cada condición experimental determinada, cada 4 milisegundos, obteniendo un total de 350 mediciones por paciente en cada uno de los 19 electrodos muestreados (dispuestos en ubicaciones establecidas por el sistema 10-20).

Sistema de referencia 10-20

El registro del EEG se realiza con distintos tipos de electrodos distribuidos según un mapa normalizado por la Federación de Electroencefalografia y Neurofisiologia Clinica, en base al sistema de coordenadas 10-20, que determina las coordenadas de cada punto de registro mediante porcentajes (10-20) de la distancia existente entre puntos del cráneo bien definidos. Esto permite el estudio de la misma región independientemente de las dimensiones y de la forma del cráneo.

En el primer Congreso Internacional de EEG, que se celebró en Londres en 1947, se reconoció que era necesario un método estándar de la colocación de los electrodos utilizados en el electroencefalograma (EEG). Varias discusiones entorno al tema dieron lugar a la definición del sistema de electrodos 10-20 (Hasper, 1958). Desde entonces, este sistema de electrodos se ha convertido en el estándar para los clínicos y para el estudio de potenciales relacionados con eventos (potenciales evocados) en entornos no clínicos.

Con el fin de determinar bajo este sistema las coordenadas de los 19 electrodos utilizados en el experimento anteriormente descrito, se utilizaron las coordenadas de todos los lugares del electrodo en una superficie de la cabeza \(real\), sobre la base de las distancias a lo largo de la superficie (triangular) de la cabeza. La superficie de la cabeza utilizada fue construida a partir de la resonancia magnética canónica que se incluye en el paquete \(SPM2\), y las ubicaciones se expresan en coordenadas \(MNI\) (standard Montreal Neurological Institute) (Oostenveld & Praamstra, 2001)

En la figura se presenta la gráfica de las coordenadas de los electrodos bajo el sistema 10-20:

Análisis de varianza a dos vías de clasificación

La energía liberada por el cerebro al ser estimulado con palabras de distinta carga emocional es diferente si el paciente es hombre o mujer. En el caso de la localización del electrodo \(C3\), los hombres presentan diferencias en algunos tiempos específicos mientras que las mujeres no. En esta sección se propone realizar una análisis de varianza con dos factores fijos que corresponden al tipo de estímulo al cual es sometido el paciente y al género del mismo.

En este caso el proceso de suavizamiento se realizó con 40 funciones base. La curva media de todas las observaciones realizadas en la localización del electrodo \(C3\) establece que el pico máximo de liberación de energía se presenta alrededor del milisegundo 400, es decir 200 milisegundos después de la aparición del estímulo.

Las funciones medias estimadas, sugieren que los hombres estimulados con palabras de carga positiva emiten menos energía y las mujeres estimuladas con palabras de carga negativa emiten mayor energía.

Al aplicar la prueba \(F\) funcional propuesta por Shen y Faraway (2004), se concluye con un nivel de significancia del 5% que no existen diferencias significativas en la respuesta bioelectrica del cerebro de mujeres y hombres cuando son estimulados con palabras de distinta carga emocional (positivos vs negativos).

#################################################
# Código para hacer  Regresión Funcional con
# los datos de un lectrodo determinado, con los  
# pacientes de ambos genero determinado, los 2 estimulos 
#(palabras positivas y negativas),para un total de
# 60 filas) y 350 mediciones en el tiempo.
##################################################

rm(list=ls())
library(zoo)
library(fda)

datos<-read.table("C:/Users/jpari/Dropbox/FDA/BrainDT.txt", header=TRUE, dec = ",")


levels(datos$esite)=c("C3","C4","CZ","F3","F4","F7","F8","FP1","FP2","FZ","O1","O2","P3","P4","PZ","T3","T4","T5","T6","VEOG")

#################################################
# Creación de la base, determinando el sexo de los
# pacientes  y el electrodo que se quiere estudiar
##################################################

BASEPOSI<-(datos$esite=="C3")&(datos$wordtype=="1")&(datos$subject !="1")

BASENEGA<-(datos$esite=="C3")&(datos$wordtype=="2")&(datos$subject !="1")
BASEPOSI = as.matrix(datos[BASEPOSI,4:353])
BASENEGA = as.matrix(datos[BASENEGA,4:353])

BASETOTAL =(datos$esite=="C3")&((datos$wordtype=="1")|(datos$wordtype=="2"))&(datos$subject !="1")
BASETOTAL =as.matrix(datos[BASETOTAL,4:353])
BASETOTALT = t(BASETOTAL)
t = dim(BASETOTALT)[1]
n = dim(BASETOTALT)[2]

BASEPOSIT = t(BASEPOSI)
dim(BASEPOSIT)

## [1] 350  46

BASENEGAT = t(BASENEGA)
dim(BASENEGAT)

## [1] 350  46

#################################################
# Creación de una base B spline de orden 4 y k = 40.
##################################################

time<- c( seq(4,1400,4))
rangeval = c(4,1400)
nbasis = 40
norder = 4
basis6 = create.bspline.basis(rangeval, nbasis, norder)

col <-matrix(2, nrow=1, ncol=n)

par(mfrow=c(1,2))
fdhposi = Data2fd(argvals = time, y = BASEPOSIT, basisobj = basis6)
plot(fdhposi,col=1, xlab="Tiempo", ylab="Voltaje", lty=1,pch = 56, main = "Estimulo positivo" )

## [1] "done"

fdhnega = Data2fd(time,BASENEGAT,  basis6)
plot(fdhnega, xlab="Tiempo", ylab="Voltaje", lty=1,pch = 56, col = 4, main = "Estimulo negativo" )

## [1] "done"

fdh = Data2fd( time,BASETOTALT, basis6)
plot(fdh,col=rep(c(1,4), n), xlab="Tiempo(milisegundos)", ylab="Voltaje(Voltios)", lty=1,pch = 56)

## [1] "done"

#################################################
#Determinación y grafica de la media y la varianza.
##################################################

meanfdh <- mean.fd(fdh)
varfdh <-var.fd(fdh)
stdvfdh  <- stddev.fd(fdh)
par(mfrow=c(1,2))

plot(meanfdh,               main="Media", xlab="Tiempo", ylab="Voltaje" )

## [1] "done"

plot(stdvfdh, main="Desviación estándar", xlab="Tiempo", ylab="Voltaje")

## [1] "done"

plot(stdvfdh^2, main="Varianza", xlab="Tiempo", ylab="Voltaje")

## [1] "done"

windows()
plot(fdh,col=8, xlab="Tiempo (milisegundos)", ylab="Voltaje (Voltios)", lty=1)

## [1] "done"

lines(meanfdh,col=2,lwd=2)
lines(stdvfdh, ylim=c(0,10), main="Desviación estándar", xlab="Tiempo (milisegundos)", ylab="Voltaje (Voltios)",lty=2, col=3, lwd=2)

#################################################
# Modelo lineal con variable respuesta la curva de
# la energia liberada por el cerebro y  variable 
# explicativa  el tipo de estímulo al cual el paciente 
# es sometido.
##################################################

matrizh <- eval.basis(time, basis6)
y2cMap <- solve(crossprod(matrizh)) %*% t(matrizh)


#################################################
#Nombre de los estímulos
##################################################


estimulos <- c("Efecto global", "Positivo  ", "Negativo","Hombre", "Mujer")

#################################################
#Indices para los estimulos y los generos
##################################################

posiindex <- c(1,   3,  5,  7,  9,  11, 13, 15, 17, 19, 21,
               23,  25, 27, 29, 31, 33, 35,37,39,41,43,45,47,49,51,53,55,57)
negindex <- c(2,    4,  6,  8,  10, 12, 14, 16, 18, 20,
              22,   24, 26, 28, 30, 32, 34, 36,38,40,42,44,46,48,50,52,54,56,58)

hombres = seq(1,58,1)
mujeres = seq(59,92,1)
#################################################
#Matriz diseño
##################################################

zmat <- matrix(0,92,5)
zmat[        ,1] <- 1
zmat[posiindex,2] <- 1
zmat[negindex,3] <- 1
zmat[hombres,4] <- 1
zmat[mujeres,5] <- 1

#################################################
#Generación de la fila, para garantizar la ortogonalidad
#y estimabilidad de los parámetros referentes a estímulos.
#Fila (0,1,1)
##################################################

z93    <- matrix(1,1,5)
z93[1] <- 0
z93[4] <- 0
z93[5] <- 0
zmat   <- rbind(zmat, z93)

#################################################
#Estimación de la fila generada a través de las 
#funciones base.
#Fila (0,1,1).
##################################################


coef   <- fdh$coefs
coef93 <- cbind(coef,matrix(0,nbasis,1))
fdh$coefs <- coef93


#################################################
#Generación de las filas, para garantizar la ortogonalidad
#y estimabilidad de los parámetros referentes a genero.
#Fila (0,1,1)
##################################################

z94    <- matrix(1,1,5)
z94[1] <- 0
z94[2] <- 0
z94[3] <- 0
zmat   <- rbind(zmat, z94)

#################################################
#Estimación de la fila generada a través de las 
#funciones base.
#Fila (0,1,1).
##################################################


coef   <- fdh$coefs
coef94 <- cbind(coef,matrix(0,nbasis,1))
fdh$coefs <- coef94


#################################################
# Estimación de los parametros funcionales, a través 
# de 20 funciones base y de orden 4
##################################################

p <- 5
xfdlist <- vector("list",p)
for (j in 1:p) xfdlist[[j]] <- zmat[,j]


nbetabasis <- 20
betabasis = create.bspline.basis(rangeval, nbetabasis, norder)

betafd    <- fd(matrix(0,nbetabasis,1), betabasis)
estimate  <- T
lambda    <- 0
betafdPar <- fdPar(betafd)

betalist <- vector("list",p)
for (j in 1:p) betalist[[j]] <- betafdPar

fRegressList <- fRegress(fdh, xfdlist, betalist)

#################################################
# Gráfica de parámetros funcionales (betas) estimados
##################################################


betaestlist <- fRegressList$betaestlist
par(mfrow=c(1,5))
for (j in 1:p) 
{
  betaestParfdj <- betaestlist[[j]]
  plot(betaestParfdj$fd, xlab="Tiempo (milisegundos)", ylab="Voltaje (Voltios)",
       main=estimulos[j])
}

#################################################
# Gráfica de funciones predichas
##################################################

yfdpar <- fRegressList$yfdpar

yhatfdobj <- fRegressList$yhatfdobj
plot(yhatfdobj)

## [1] "done"

#################################################
# Cálculo de la matriz de residuales y obtención de
# la matriz de varianza y covarianza
##################################################


yhatmat  <- eval.fd(time,yhatfdobj)      # Valores estimados
matplot(yhatmat[,1:92], col=2, type="l")
ymat     <- eval.fd(time, fdh)              # Valores observados 
matplot(time,ymat[,1:92], type="l")
alfa <- betaestlist[[1]]                              
alpha    <-  eval.fd(time,alfa$fd)          # Constantes (Media) estimada
plot(alfa$fd, xlab="Tiempo", ylab="Voltaje")

## [1] "done"

Residuales <- ymat[,1:92] - yhatmat[,1:92]  # Los residuales son iguales a observados - estimados
matplot(Residuales, type = "l")

##########################################################################
# Gráfica de contorno de la varianza de los residuales del modelo completo
##########################################################################


SigmaE   <- var(t(Residuales))
contour(SigmaE, xlab="Milisegundos", ylab="Milisegundos", cex=1.2)
dim(SigmaE)

## [1] 350 350

######################################################################
# Gráfica de dsviación estandar de los residuales del modelo completo
######################################################################
cexval = 1.2
par(mfrow=c(1,1), mar=c(5,5,4,2)+cexval+2, pty="m")

stddevE <- sqrt(diag(SigmaE))
plot(time, stddevE, type="l", cex=1.2,
     xlab="Milisegundos", ylab="Desviacion estandar (Voltios)")

##########################################################################
# Estimación del modelo reducido 1: Solo los estímulos
##########################################################################
#################################################
#Nombre de los estímulos
##################################################


estimulos <- c("Efecto global", "Positivo  ", "Negativo")

#################################################
#Indices para los estimulos 
##################################################

posiindex1 <- seq(1,92,2)
negindex1 <- seq(2,92,2)

#################################################
#Matriz diseño
##################################################
q1 = 3
fdh1 = Data2fd(time,BASETOTALT, basis6)
zmat1 <- matrix(0,92,3)
zmat1[        ,1] <- 1
zmat1[posiindex1,2] <- 1
zmat1[negindex1,3] <- 1

#################################################
#Generación de la fila, para garantizar la ortogonalidad
#y estimabilidad de los parámetros.
#Fila (0,1,1,1)
##################################################

z931    <- matrix(1,1,3)
z931[1] <- 0
zmat1   <- rbind(zmat1, z931)

#################################################
#Estimación de la fila generada a través de las 
#funciones base.
#Fila (0,1,1,1).
##################################################


coef   <- fdh1$coefs
coef931 <- cbind(coef,matrix(0,nbasis,1))
fdh1$coefs <- coef931

#################################################
# Estimación de los parametros funcionales, a través 
# de 20 funciones base y de orden 4
##################################################

xfdlist1 <- vector("list",q1)
for (j in 1:q1) xfdlist1[[j]] <- zmat1[,j]


nbetabasis <- 20
betabasis = create.bspline.basis(rangeval, nbetabasis, norder)

betafd    <- fd(matrix(0,nbetabasis,1), betabasis)
estimate  <- T
lambda    <- 0
betafdPar <- fdPar(betafd)

betalist1 <- vector("list",q1)
for (j in 1:q1) betalist1[[j]] <- betafdPar

fRegressList1 <- fRegress(fdh1, xfdlist1, betalist1)

#################################################
# Gráfica de parámetros funcionales (betas) estimados
##################################################


betaestlist1 <- fRegressList1$betaestlist
par(mfrow=c(3,1))
for (j in 1:q1) 
{
  betaestParfdj <- betaestlist1[[j]]
  plot(betaestParfdj$fd, xlab="Tiempo (milisegundos)", ylab="Voltaje (Voltios)",
       main=estimulos[j])
}

#################################################
# Gráfica de funciones predichas
##################################################

yfdpar1 <- fRegressList1$yfdpar

yhatfdobj1 <- fRegressList1$yhatfdobj
plot(yhatfdobj1)

## [1] "done"

#################################################
# Cálculo de la matriz de residuales y obtención de
# la matriz de varianza y covarianza del modelo reducido 1: Solo estimulos
##################################################


yhatmat1  <- eval.fd(time,yhatfdobj1)      # Valores estimados
matplot(yhatmat1[,1:n], col=2, type="l")
ymat1     <- eval.fd(time, fdh1)              # Valores observados 
matplot(time,ymat1[,1:n], type="l")

alfa <- betaestlist1[[1]]                              
alpha    <-  eval.fd(time,alfa$fd)          # Constantes (Media) estimada
plot(alfa$fd, xlab="Tiempo", ylab="Voltaje")

## [1] "done"

Residuales1 <- ymat1[,1:n] - yhatmat1[,1:n]  # Los residuales son iguales a observados - estimados
matplot(Residuales1, type = "l")


SigmaRed1E   <- var(t(Residuales1))
contour(SigmaRed1E, xlab="Milisegundos", ylab="Milisegundos", cex=1.2)

dim(SigmaRed1E)

## [1] 350 350

######################################################################
# Gráfica de dsviación estandar de los residuales del modelo reducido
######################################################################
cexval = 1.2
par(mfrow=c(1,1), mar=c(5,5,4,2)+cexval+2, pty="m")
stddevRed1E <- sqrt(diag(SigmaRed1E))
plot(time, stddevRed1E, type="l", cex=1.2,
     xlab="Milisegundos", ylab="Desviacion estandar (Voltios)")

######################################################################
# Prueba F de Fareway para comparar modelo lleno menos modelo reducido
######################################################################

Num =( sum(diag(SigmaRed1E-SigmaE))/(p-q1))
Den = (sum(diag(SigmaE))/(n-p))
F1 = Num/Den
factor = (((sum(diag(SigmaE)))^2)/ (sum(diag(SigmaE^2))))
f1 = factor*(p-q1)
f2 = factor*(n-p)
Fteo1<-qf(0.95,f1,f2)



##########################################################################
# Estimación del modelo reducido 2: Solo los géneros
##########################################################################
#################################################
#Nombre de los estímulos
##################################################
q2 = 3
fdh2 = Data2fd(time, BASETOTALT, basis6)
estimulos <- c("Efecto global", "Hombres", "Mujeres")

#################################################
#Indices para los estimulos 
##################################################

Hombres1 <- seq(1,58,1)
Mujeres1 <- seq(59,92,1)

#################################################
#Matriz diseño
##################################################

zmat2 <- matrix(0,92,3)
zmat2[        ,1] <- 1
zmat2[Hombres1,2] <- 1
zmat2[Mujeres1,3] <- 1

#################################################
#Generación de la fila, para garantizar la ortogonalidad
#y estimabilidad de los parámetros.
#Fila (0,1,1,1)
##################################################

z932   <- matrix(1,1,3)
z932[1] <- 0
zmat2   <- rbind(zmat2, z932)

#################################################
#Estimación de la fila generada a través de las 
#funciones base.
#Fila (0,1,1,1).
##################################################


coef   <- fdh2$coefs
coef932 <- cbind(coef,matrix(0,nbasis,1))
fdh2$coefs <- coef932

#################################################
# Estimación de los parametros funcionales, a través 
# de 20 funciones base y de orden 4
##################################################

xfdlist2 <- vector("list",q2)
for (j in 1:q2) xfdlist2[[j]] <- zmat2[,j]


nbetabasis <- 20
betabasis = create.bspline.basis(rangeval, nbetabasis, norder)

betafd    <- fd(matrix(0,nbetabasis,1), betabasis)
estimate  <- T
lambda    <- 0
betafdPar <- fdPar(betafd)

betalist2 <- vector("list",q2)
for (j in 1:q2) betalist2[[j]] <- betafdPar

fRegressList2 <- fRegress(fdh2, xfdlist2, betalist2)


#################################################
# Gráfica de parámetros funcionales (betas) estimados
##################################################


betaestlist2 <- fRegressList2$betaestlist
par(mfrow=c(3,1))
for (j in 1:q2) 
{
  betaestParfdj <- betaestlist2[[j]]
  #plot(betaestParfdj$fd, xlab="Tiempo (milisegundos)", ylab="Voltaje (Voltios)", main=estimulos[j])
}


#################################################
# Gráfica de funciones predichas
##################################################

yfdpar2 <- fRegressList2$yfdpar

yhatfdobj2 <- fRegressList2$yhatfdobj
#plot(yhatfdobj2)




#################################################
# Cálculo de la matriz de residuales y obtención de
# la matriz de varianza y covarianza del modelo reducido 1: Solo estimulos
##################################################


yhatmat2  <- eval.fd(time,yhatfdobj2)      # Valores estimados
#matplot(yhatmat2[,1:n], col=2, type="l")
ymat2     <- eval.fd(time, fdh)              # Valores observados 
#matplot(time,ymat2[,1:n], type="l")
alfa <- betaestlist2[[1]]                              
alpha    <-  eval.fd(time,alfa$fd)          # Constantes (Media) estimada
#plot(alfa$fd, xlab="Tiempo", ylab="Voltaje")
Residuales2 <- ymat2[,1:n] - yhatmat2[,1:n]  # Los residuales son iguales a observados - estimados
#matplot(Residuales2, type = "l")


SigmaRed2E   <- var(t(Residuales2))
#contour(SigmaRed2E, xlab="Milisegundos", ylab="Milisegundos", cex=1.2)
dim(SigmaRed2E)

## [1] 350 350

######################################################################
# Gráfica de dsviación estandar de los residuales del modelo reducido
######################################################################
cexval = 1.2
par(mfrow=c(1,1), mar=c(5,5,4,2)+cexval+2, pty="m")
stddevRed2E <- sqrt(diag(SigmaRed2E))
plot(time, stddevRed2E, type="l", cex=1.2,
     xlab="Milisegundos", ylab="Desviacion estandar (Voltios)")

######################################################################
# Prueba F de Fareway para comparar modelo lleno menos modelo reducido
######################################################################
Num =( sum(diag(SigmaRed2E-SigmaE))/(p-q2))
Den = (sum(diag(SigmaE))/(n-p))
F2 = Num/Den
factor = (((sum(diag(SigmaE)))^2)/ (sum(diag(SigmaE^2))))
f12 = factor*(p-q2)
f22 = factor*(n-p)
Fteo2<-qf(0.95,f12,f22)