Analisis_baza_006.Rmd

---
title: "Analisis_baza_006"
output: word_document
---
# Cargar paquetes
```{r, include=FALSE}
#install.packages("lubridate")
library(lubridate) # Para las fechas
```

# DATOS INSTAR

## Importar datos INSTAR 

- La tabla de datos se construye en excel y se guarda como .csv
- Se abre con procesador de texto, y se comprueba que la fecha esta en formato mm-dd-yy
- Se guarda como .txt
 Las columnas que constituyen la tabla son: fecha, huevo, L1, L2, crisalida, radiacion, tmax, tmin, tmed y exergia.

```{r, echo=FALSE}
INSTAR<-read.table("instar_bz006.txt", header=T, sep=",", dec=".")
INSTAR$fecha <-as.Date(INSTAR$fecha, format="%m-%d-%y")
```

## Anadir biociclos al data.frame
```{r, echo=FALSE}
INSTAR$yday <- yday(INSTAR$fecha) # anadimos el dia del anio con el paquete lubridate
INSTAR$year <- year(INSTAR$fecha) # anadimos el anio con el paquete lubridate
INSTAR$biociclo <- NULL
#write.csv(INSTAR, file = "INSTAR_baza006.csv", row.names=FALSE, na="")
```

Aparece dos chunks para este proceso, uno para Baza y otro para Sierra Nevada, ya que tienen fechas de inicio diferentes. 

```{r, echo=FALSE}
#Baza
subset_1993 <- subset(INSTAR, INSTAR$year==1993)
for(i in 1:length(subset_1993$yday)){
  ifelse(subset_1993$yday[i]<=79, 
       subset_1993$biociclo[i] <- 1, subset_1993$biociclo[i] <- 2)
}

subset_1994 <- subset(INSTAR, INSTAR$year==1994)
for(i in 1:length(subset_1994$yday)){
  ifelse(subset_1994$yday[i]<=79, 
       subset_1994$biociclo[i] <- 2, subset_1994$biociclo[i] <- 3)
}

subset_1995 <- subset(INSTAR, INSTAR$year==1995)
for(i in 1:length(subset_1995$yday)){
  ifelse(subset_1995$yday[i]<=79, 
       subset_1995$biociclo[i] <- 3, subset_1995$biociclo[i] <- 4)
}

subset_1996 <- subset(INSTAR, INSTAR$year==1996)
for(i in 1:length(subset_1996$yday)){
  ifelse(subset_1996$yday[i]<=79, 
       subset_1996$biociclo[i] <- 4, subset_1996$biociclo[i] <- 5)
}

subset_1997 <- subset(INSTAR, INSTAR$year==1997)
for(i in 1:length(subset_1997$yday)){
  ifelse(subset_1997$yday[i]<=79, 
       subset_1997$biociclo[i] <- 5, subset_1997$biociclo[i] <- 6)
}

subset_1998 <- subset(INSTAR, INSTAR$year==1998)
for(i in 1:length(subset_1998$yday)){
  ifelse(subset_1998$yday[i]<=79, 
       subset_1998$biociclo[i] <- 6, subset_1998$biociclo[i] <- 7)
}

subset_1999 <- subset(INSTAR, INSTAR$year==1999)
for(i in 1:length(subset_1999$yday)){
  ifelse(subset_1999$yday[i]<=79, 
       subset_1999$biociclo[i] <- 7, subset_1999$biociclo[i] <- 8)
}

subset_2000 <- subset(INSTAR, INSTAR$year==2000)
for(i in 1:length(subset_2000$yday)){
  ifelse(subset_2000$yday[i]<=79, 
       subset_2000$biociclo[i] <- 8, subset_2000$biociclo[i] <- 9)
}

subset_2001 <- subset(INSTAR, INSTAR$year==2001)
for(i in 1:length(subset_2001$yday)){
  ifelse(subset_2001$yday[i]<=79, 
       subset_2001$biociclo[i] <- 9, subset_2001$biociclo[i] <- 10)
}

subset_2002 <- subset(INSTAR, INSTAR$year==2002)
for(i in 1:length(subset_2002$yday)){
  ifelse(subset_2002$yday[i]<=79, 
       subset_2002$biociclo[i] <- 10, subset_2002$biociclo[i] <- 11)
}

subset_2003 <- subset(INSTAR, INSTAR$year==2003)
for(i in 1:length(subset_2003$yday)){
  ifelse(subset_2003$yday[i]<=79, 
       subset_2003$biociclo[i] <- 11, subset_2003$biociclo[i] <- 12)
}

subset_2004 <- subset(INSTAR, INSTAR$year==2004)
for(i in 1:length(subset_2004$yday)){
  ifelse(subset_2004$yday[i]<=79, 
       subset_2004$biociclo[i] <- 12, subset_2004$biociclo[i] <- 13)
}
subset_2005 <- subset(INSTAR, INSTAR$year==2005)
for(i in 1:length(subset_2005$yday)){
  ifelse(subset_2005$yday[i]<=79, 
       subset_2005$biociclo[i] <- 13, subset_2005$biociclo[i] <- 14)
}

subset_2006 <- subset(INSTAR, INSTAR$year==2006)
for(i in 1:length(subset_2006$yday)){
  ifelse(subset_2006$yday[i]<=79, 
       subset_2006$biociclo[i] <- 14, subset_2006$biociclo[i] <- 15)
}

subset_2007 <- subset(INSTAR, INSTAR$year==2007)
for(i in 1:length(subset_2007$yday)){
  ifelse(subset_2007$yday[i]<=79, 
       subset_2007$biociclo[i] <- 15, subset_2007$biociclo[i] <- 16)
}

subset_2008 <- subset(INSTAR, INSTAR$year==2008)
for(i in 1:length(subset_2008$yday)){
  ifelse(subset_2008$yday[i]<=79, 
       subset_2008$biociclo[i] <- 16, subset_2008$biociclo[i] <- 17)
}

subset_2009 <- subset(INSTAR, INSTAR$year==2009)
for(i in 1:length(subset_2009$yday)){
  ifelse(subset_2009$yday[i]<=79, 
       subset_2009$biociclo[i] <- 17, subset_2009$biociclo[i] <- 18)
}

subset_2010 <- subset(INSTAR, INSTAR$year==2010)
for(i in 1:length(subset_2010$yday)){
  ifelse(subset_2010$yday[i]<=79, 
       subset_2010$biociclo[i] <- 18, subset_2010$biociclo[i] <- 19)
}

subset_2011 <- subset(INSTAR, INSTAR$year==2011)
for(i in 1:length(subset_2011$yday)){
  ifelse(subset_2011$yday[i]<=79, 
       subset_2011$biociclo[i] <- 19, subset_2011$biociclo[i] <- 20)
}

subset_2012 <- subset(INSTAR, INSTAR$year==2012)
for(i in 1:length(subset_2012$yday)){
  ifelse(subset_2012$yday[i]<=79, 
       subset_2012$biociclo[i] <- 20, subset_2012$biociclo[i] <- 21)
}

subset_2013 <- subset(INSTAR, INSTAR$year==2013)
for(i in 1:length(subset_2013$yday)){
  ifelse(subset_2013$yday[i]<=79, 
       subset_2013$biociclo[i] <- 21, subset_2013$biociclo[i] <- 22)
}

subset_2014 <- subset(INSTAR, INSTAR$year==2014)
for(i in 1:length(subset_2014$yday)){
  ifelse(subset_2014$yday[i]<=79, 
       subset_2014$biociclo[i] <- 22, subset_2014$biociclo[i] <- 23)
}

INSTAR<- rbind(subset_1994, subset_1995, subset_1996, subset_1997, subset_1998, subset_1999, subset_2000, subset_2001, subset_2002, subset_2003, subset_2004, subset_2005, subset_2006, subset_2007, subset_2008, subset_2009, subset_2010, subset_2011, subset_2012, subset_2013, subset_2014) #Se quitan el primer y ultimo biociclo porque no estan completos

#QUITAR BIOCICLO 2 y 23(incompletos)
INSTAR<-INSTAR[!(INSTAR$biociclo=="2"),]
INSTAR<-INSTAR[!(INSTAR$biociclo=="23"),]

```

## Agregacion INSTAR

Media del vigor y sumas (l1, l2 y l1+l2) por biociclo

```{r, echo=FALSE}
sum_l1<- aggregate(x = INSTAR$L1, by=list(biociclo=INSTAR$biociclo), FUN=sum, na.rm=TRUE)
names(sum_l1)[2] <- "sum_l1"

sum_l2<- aggregate(x = INSTAR$L2, by=list(biociclo=INSTAR$biociclo), FUN=sum, na.rm=TRUE)
names(sum_l2)[2] <- "sum_l2"

larvas<-cbind(sum_l1, sum_l2)
larvas<-larvas[-c(3)]
larvas$sum_l1l2=larvas[,2]+larvas[,3]

sum_huevo<- aggregate(x = INSTAR$huevo, by=list(biociclo=INSTAR$biociclo), FUN=sum, na.rm=TRUE)
names(sum_huevo)[2] <- "sum_huevo"

sum_crisalida<- aggregate(x = INSTAR$crisalida, by=list(biociclo=INSTAR$biociclo), FUN=sum, na.rm=TRUE)
names(sum_crisalida)[2] <- "sum_crisalida"

avg_exergia<-aggregate(INSTAR$exergia, by=list(biociclo=INSTAR$biociclo), FUN=mean, na.rm=TRUE)
names(avg_exergia)[2] <- "avg_exergia"


agreg_INSTAR <- cbind(larvas,avg_exergia, sum_crisalida, sum_huevo)
agreg_INSTAR_biociclo <- agreg_INSTAR[ -c(5, 7, 9) ]

#write.csv(agreg_INSTAR_biociclo,"agreg_INSTAR_bz008_biociclo.csv",row.names=FALSE, na="")
```

# DATOS COPLAS

## Importar datos COPLAS 

```{r}
COPLAS<-read.csv("zona_baza.csv", header=TRUE, sep=",")
names(COPLAS)[1] <- "RODAL"
```

Activar el rodal necesario para el analisis

```{r, echo=FALSE}
rodal<-subset(COPLAS, RODAL == "GR000006")
#rodal<-subset(COPLAS, RODAL == "GR023008")
#rodal<-subset(COPLAS, RODAL == "GR023012")
#rodal<-subset(COPLAS, RODAL == "GR023015")
#rodal<-subset(COPLAS, RODAL == "GR023017")
#rodal<-subset(COPLAS, RODAL == "GR023021")
```

```{r, echo=FALSE}
#Eliminar primer y ultimo anio porque en Instar no teneos datos de ese anio completo. 
rodal<-rodal[-c(1,2, 23),]
```

## Definir biociclo COPLAS 
```{r, echo=FALSE}
rodal$biociclo <- 3:22 #Sierra de Baza
```

##Union de tablas INSTAR y COPLAS
En esta tabla aparecen los datos de Instar agregados, los Ln de sumL1, sumL2 y sumL1+sumL2, el anio y el grado de infestacion de COPLAS. 
```{r, echo=FALSE}
COPLAS_INSTAR<-merge.data.frame(agreg_INSTAR_biociclo, rodal, by.x = "biociclo")
COPLAS_INSTAR<-COPLAS_INSTAR[-c(8)] #Eliminar la columna con el nombre del rodal
```


# VALIDACION EXTERNA

## Correlacion [L1, L2, L1+L2 vs. COPLAS]

Ya que es dificil explicar que una variable es dependiente de otra. La mejor forma de validar el modelo (con datos reales) es haciendo una correlacion. Con este analisis sabremos si existe relacion o no entre la variable Instar y COPLAS. Sabremos tambien la intensidad de esta relacion y si es significativa.
Como nuestros datos no siguen una distribucion normal, haremos una correlacion con datos no parametricos: 'Spearman'

Funcion en R (pedir a Antonio este codigo) que plotea la relacion entre ellos con un circulo (cuanto mas relacionados, mas grande el circulo) y si es negativa o positiva. Da tambien el p valor de la correlacion. *"Puedo tener fuerte relacion pero el coef no es diferente de 0, a lo mejor con mas datos si."*" 

Deberiamos esperar que haya mas relacion entre coplas y exergia, luego coplas y larvas. Cuando hay relacion negativa es que es inversa. Y como no sabemos quien es x e y, quiere decir que cuando una aumenta, otra disminuye, peor no sabemos quien aumente y quien disminuye. 

```{r, echo=FALSE}
library(corrplot)
```

```{r, echo=FALSE}
#Preparamos una tabla con las variables que se van a correlacionar
tabla_cor<-COPLAS_INSTAR[-c(1, 6, 7, 8)]
names(tabla_cor)[5]<-"Grado infestacion"
#Hacemos la correlacion
correlacion<-cor(tabla_cor, method = "spearman")
#Calculo significancia
cor.mtest <- function(mat, conf.level = 0.95){
  mat <- as.matrix(mat)
    n <- ncol(mat)
    p.mat <- lowCI.mat <- uppCI.mat <- matrix(NA, n, n)
    diag(p.mat) <- 0
    diag(lowCI.mat) <- diag(uppCI.mat) <- 1
    for(i in 1:(n-1)){
        for(j in (i+1):n){
            tmp <- cor.test(mat[,i], mat[,j], conf.level = conf.level)
            p.mat[i,j] <- p.mat[j,i] <- tmp$p.value
            lowCI.mat[i,j] <- lowCI.mat[j,i] <- tmp$conf.int[1]
            uppCI.mat[i,j] <- uppCI.mat[j,i] <- tmp$conf.int[2]
        }
    }
    return(list(p.mat, lowCI.mat, uppCI.mat))
}
res1 <- cor.mtest(tabla_cor,0.95)
```


```{r}
#Plot
corrplot(correlacion, p.mat = res1[[1]], sig.level=0.05, method = "number", type = "upper", tl.col = "black", number.cex = 1.5, insig = "p-value")
```

- La correlacion de exergia y coplas resulta ser menos fiable de los esperado. Posiblemente porque su modelizacion es simple.

- La correlacion frente a COPLAS de L1+L2 y L1 son parecidas. Esto quiere decir que L2 no aporta mucho (o nada) cuando se agregan. Esto probablemente sea porque la agregacion de las larvas indica la carga defoliadora en el sistema y al hacer esta agregacion (pobre) se le da la misma capacidad defoliadora a L1 que a L2, cuando en realidad L2 tiene mas. 

# VALIDACION INTERNA

##ACP
Autocorrelation plot. 
Nos ayudara a saber si existen ciclos dentro de toda la serie temporal. 
ESTE ANALISIS SOLO SE HACE CON LOS DATOS DE BAZA

"*Those plots are showing you the correlation of the series with itself, lagged by x time unitscorrelation of the series with itself, lagged by x time units. So imagine taking your time series of length TT, copying it, and deleting the first observation of copy#1 and the last observation of copy#2. Now you have two series of length T for which you calculate a correlation coefficient. This is the value of of the vertical axis at x=1 in your plots. It represents the correlation of the series lagged by one time unit. You go on and do this for all possible time lags xx and this defines the plot*.

*The answer to your question of what is needed to report a pattern is dependent on what pattern you would like to report. But quantitatively speaking, you have exactly what I just described: the correlation coefficient at different lags of the series. You can extract these numerical values by issuing the command  acf(x.ts,100)$acf*.

*In terms of what lag to use, this is again a matter of context. It is often the case that there will be specific lags of interest. Say, for example, you may believe the fish species migrates to and from an area every ~30 days. This may lead you to hypothesize a correlation in the time series at lags of 30. In this case, you would have support for your hypothesis*"

"*The blue lines give the values beyond which the autocorrelations are (statistically) significantly different from zero*" 

Hay que definir el LAG. lag.max= Lo hemos definido como el numero de biociclos.

**CORRELOGRAMAS INSTAR**
```{r, echo=FALSE}
exergia_ts_y<-ts(agreg_INSTAR_biociclo$avg_exergia, frequency = 1, start = c(1994,84), end = c(2014, 79))
exergia_acf_y<-acf(exergia_ts_y, type=c("correlation"), plot=TRUE, lag.max = 19)

L1_ts_y<-ts(agreg_INSTAR_biociclo$sum_l1, frequency = 1, start = c(1994,84), end = c(2014, 79))
L1_acf_y<-acf(L1_ts_y, type=c("correlation"), plot=TRUE, lag.max = 19)

L2_ts_y<-ts(agreg_INSTAR_biociclo$sum_l2, frequency = 1, start = c(1994,84), end = c(2014, 79))
L2_acf_y<-acf(L2_ts_y, type=c("correlation"), plot=TRUE, lag.max = 19)

huevo_ts_y<-ts(agreg_INSTAR_biociclo$sum_huevo, frequency = 1, start = c(1994,84), end = c(2014, 79))
huevo_acf_y<-acf(huevo_ts_y, type=c("correlation"), plot=TRUE, lag.max = 19)

crisalida_ts_y<-ts(agreg_INSTAR_biociclo$sum_crisalida, frequency = 1, start = c(1994,84), end = c(2014, 79))
crisalida_acf_y<-acf(crisalida_ts_y, type=c("correlation"), plot=TRUE, lag.max = 19)
```


**CORRELOGRAMA COPLAS**
```{r, echo=FALSE}
#Mirar ACF de COPLAS
coplas_ts<-ts(COPLAS$GRADO.REVISADO, frequency = 1)
coplas_acf<-acf(coplas_ts, type=c("correlation"), plot=TRUE, lag.max = 19)
```

Los Lags 1, 5, 6 y 17 son significativos. ??
La de 5 y 6 podria referirse al ciclo del que hablan en la literatura. 

##Analisis de tendencia

```{r}
ind_feno<-read.csv("indicadores_feno_bz006.csv", header = TRUE, sep = ",", dec = ".")
```

```{r}
library(Kendall)
```

Se eligen indicadores de cada variable para calcular su tendencia.

En las variables huevo, L1, L2 y crisalida:

- Fecha del maximo. 

- Fecha del primer dia distinto a cero. 

- Fecha del ultimo dia distinto a cero. 


En la variable exergia:

- Fecha del maximo.

- Fecha del minimo. 

H0= no ha tendencia. Si p<0.05 entonces si hay tendencia!

**TENDENCIA EN HUEVOS**

*DIA DE PICO MAXIMO*
```{r, echo=FALSE}
max_huevo <- merge(aggregate(huevo ~ biociclo, data = INSTAR, FUN = max), INSTAR)
max_huevo <- max_huevo[ -c(2:11,13) ] #Quedarme solo con biociclo-dia
names(max_huevo)[2] <- "dia_max_huevo"
tend_max_huevo<-MannKendall(max_huevo$dia_max_huevo)
summary(tend_max_huevo)
plot(max_huevo, ylab="Dia con el maximo numero de huevos", las=1, xlab="Tiempo")
```

Los valores del tau de MannKendall van de -1 a 1, por lo tanto un tau =0.174 no indicaria tendencia. Aunque no podriamos asegurarlo ya que el dato no es significativo. Muy marcada la diferencia de los 3 primeros datos con el resto, podra deberse al efecto llamado 'calentamiento dle modelo' Tal vez, podriamos decir que el modelo en el tiempo se va 'adecuando' a las fechas reales, si es que la fecha real del pico max de huevo es sobre el dia 40 (mediados de febrero). 

*PRIMER DIA*
```{r, echo=FALSE}
tend_ini_huevo<-MannKendall(ind_feno$I_huevo)
summary(tend_ini_huevo)
plot(ind_feno$I_huevo, ylab="Dia del inicio del periodo de presencia de huevos", las=1, xlab="Tiempo")
```

Tendencia negativa (tau= -0.43) y significativa. Puede deberse al calentamiento dle modelo..

*ULTIMO DIA*
```{r, echo=FALSE}
tend_fin_huevo<-MannKendall(ind_feno$F_huevo)
summary(tend_fin_huevo)
plot(ind_feno$F_huevo, ylab="Dia del final del periodo de presencia de huevos", las=1, xlab="Tiempo")
```

Tendencia casi nula pero esto no es significativo

**TENDENCIA EN L1**

*DIA DE PICO MAXIMO*
```{r, echo=FALSE}
max_l1 <- merge(aggregate(L1 ~ biociclo, data = INSTAR, FUN = max), INSTAR, raw.names=FALSE)
max_l1 <- max_l1[ -c(2:11,13) ] #Quedarme solo con biociclo-dia
names(max_l1)[2] <- "dia_max_L1"
tend_max_l1<-MannKendall(max_l1$dia_max_L1)
summary(tend_max_l1)
plot(max_l1, ylab="Dia con el maximo numero de larvas 1", las= 1, xlab = "Tiempo")
```

Tendencia casi nula pero no es significativo. 

*PRIMER DIA*
```{r, echo=FALSE}
tend_ini_L1<-MannKendall(ind_feno$I_l1)
summary(tend_ini_L1)
plot(ind_feno$I_l1, ylab="Dia del inicio del periodo de presencia de larvas 1", las=1, xlab="Tiempo")
```

Tendencia negativa y significativa. A?os de calentamiento!

*ULTIMO DIA*
```{r, echo=FALSE}
tend_fin_L1<-MannKendall(ind_feno$F_l1)
summary(tend_fin_L1)
plot(ind_feno$F_l1, ylab="Dia del final del periodo de presencia de larvas 1", las=1, xlab="Tiempo")
```

El valor de la intensidad de la tendencia es bajo pero este dato no es significativo. 

**TENDENCIA EN L2**

*DIA DE PICO MAXIMO*
```{r, echo=FALSE}
max_l2 <- merge(aggregate(L2 ~ biociclo, data = INSTAR, FUN = max), INSTAR)
max_l2 <- max_l2[ -c(2:11,13) ] #Quedarme solo con biociclo-dia
names(max_l2)[2] <- "dia_max_L2"
tend_max_l2<-MannKendall(max_l2$dia_max_L2)
summary(tend_max_l2)
plot(max_l2, ylab="Dia con el maximo numero de larvas 2", las= 1, xlab = "Tiempo")
```

El valor de tau es cercano a 0 pero no es significativo; los datos son muy dispares.

*PRIMER DIA*
```{r, echo=FALSE}
tend_ini_L2<-MannKendall(ind_feno$I_l2)
summary(tend_ini_L2)
plot(ind_feno$I_l2, ylab="Dia del inicio del periodo de presencia de larvas 2", las=1, xlab="Tiempo")
```

Tendencia negativa y significativa debido al calentamiento del modelo. 

*ULTIMO DIA*
```{r, echo=FALSE}
tend_fin_L2<-MannKendall(ind_feno$F_l2)
summary(tend_fin_L2)
plot(ind_feno$F_l2, ylab="Dia del final inicio del periodo de presencia de larvas 2", las=1, xlab="Tiempo")
```

Datos muy dispares

**TENDENCIA EN CRISALIDA**

*DIA DE PICO MAXIMO*
```{r, echo=FALSE}
tend_max_crisalida<-MannKendall(ind_feno$max_crisalida)
summary(tend_max_crisalida)
plot(ind_feno$max_crisalida, ylab="Dia con el maximo numero de crisalidas", las=1, xlab="Tiempo")
```

Tendencia debido al calentamiento.

*PRIMER DIA*
```{r, echo=FALSE}
tend_ini_crisalida<-MannKendall(ind_feno$I_crisalida)
summary(tend_ini_crisalida)
plot(ind_feno$I_crisalida, ylab="Dia del inicio del periodo de presencia de crisalidas", las=1, xlab="Tiempo")
```

Tendencia negativa significativa. Calentamiento

*ULTIMO DIA*
```{r, echo=FALSE}
tend_fin_crisalida<-MannKendall(ind_feno$F_crisalida)
summary(tend_fin_crisalida)
plot(ind_feno$F_crisalida, ylab="Dia del final del periodo de presencia de crisalidas", las=1, xlab="Tiempo")
```

Valor de la tendencia cercano a cero, pero no significativo. Datos muy dispares. 

**TENDENCIA EN EXERGIA**

*DIA DLE PICO MAXIMO*
```{r, echo=FALSE}
max_exergia <- merge(aggregate(exergia ~ biociclo, data = INSTAR, FUN = max), INSTAR)
max_exergia <- max_exergia[ -c(2:11,13) ] #Quedarme solo con biociclo-dia
names(max_exergia)[2] <- "dia_max_exergia"
tend_max_exergia<-MannKendall(max_exergia$dia_max_exergia)
summary(tend_max_exergia)
plot(max_exergia, ylab="Dia con la maxima exergia", las= 1, xlab = "Tiempo")
```

Datos muy dispares en algunos a?os, por eso no sale significativo. 

*DIA DEL PICO MINIMO*
```{r, echo=FALSE}
min_exergia <- merge(aggregate(exergia ~ biociclo, data = INSTAR, FUN = min), INSTAR)
min_exergia <- min_exergia[ -c(2:11,13) ] #Quedarme solo con biociclo-dia
names(min_exergia)[2] <- "dia_min_exergia"
tend_min_exergia<-MannKendall(min_exergia$dia_min_exergia)
summary(tend_min_exergia)
plot(min_exergia, ylab="Dia con la minima exergia", las= 1, xlab = "Tiempo")
```

Calentamiento del modelo.