#Install and load the packages below as needed
install.packages("fda")
install.packages("car")
install.packages("MESS")
install.packages("ggplot2")

library(fda)
library(car)
library(MESS)
library(ggplot2)

setwd ("???") #set working directory and save the provided data files in it

#Part 1
########################################################################
########## Smoothing NEFA response #####################################
########################################################################
#This part of the script is adapted to smooth time series and extract dscriptors from the smoothed functions of physiological variables that increase in face to the perturbation.

cat("\014") #Clear console
rm(list=ls(all=TRUE)) #Remove the history 
dev.off() #To close all graphs

DF <- read.table("sample_cow_nefa.csv", header = TRUE,sep = ';') 
#This dataset refers to nefa response of a cow exposed to one feed restriction of 3 days (from t= 1 d until t = 4 d).
#Dataset is available for download under https://doi.org/10.5281/zenodo.5776521.

########################################################################
########## Using FDA ###################################################
########################################################################

#Smoothing coefficient and grid parameters
lambda = 1e-3 #Coefficient of roughness or penalty
range = c(min(DF$t),max(DF$t))
break.vec = DF$t[seq(1, length(DF$t),1)]#A vector with the length of the dataset especifying the knots where piecewise polinomials will be fitted.

#Basis function
basis.func = create.bspline.basis(rangeval = range, norder = 6, breaks = break.vec) #The parameter "norder" is The order of b-splines: cubic, quadratic, etc.
param.func = fdPar(fdobj = basis.func, Lfdobj = 4,lambda = lambda) #The parameter "Ldfobj" is a linear differential operator object.
smooth.nefa = smooth.basis(argvals = DF$t,DF$nefa,param.func)
nefa.fda = smooth.nefa$fd

#Estimation of response function and its first derivative
eval.days = seq(min(DF$t),max(DF$t),0.1)
nefa.estim   = eval.fd(eval.days,nefa.fda) #Estimation of nefa response
nefa.estim.1 = eval.fd(eval.days,nefa.fda,1) #Estimation of first derivative from nefa response

########################################################################
########## Extracting features from challenge ##########################
########################################################################
#Challenge duration is considered as the sum of the durations of the restriction and recovery periods.
#Determine baseline before the restriction begins
TimePpre<-1 #This is the known time when restriction started. In this example, restriction started at t = 1d

#Baseline is assumed to be constant, and determined as an average of all experimental data prior to beginning of restriction. In this case, an average of data points between t = -14d and t=-1d
Ppre<-mean((subset(DF, t> -14 & t<= -1, select=c(nefa)))$nefa) 

#Create a vector containing time and the estimation of smoothed curve corresponding to animal response to restriction
estimated_data<-data.frame(c(eval.days),c(nefa.estim))
nrow_estimated_data<-nrow(estimated_data) #Determine the number of lines of the 'estimated_data' vector.

#Create a vector containing time and the estimated baseline. This vector should have the same size as 'estimated_data' vector.
estimated_baseline1<-data.frame(rep(Ppre, nrow_estimated_data))

#Create a vector containing time and the estimation of 1st derivative from smoothed curve
der1<-data.frame(c(eval.days),c(nefa.estim.1))


#Define inflection point during  restriction
#In this example, the restriction started at t=1 d and ended at t=4 d.
c1 <- subset(der1, c.eval.days.> 1 & c.eval.days.<= 4, select=c(c.nefa.estim.1.)) #Create a vector 'c1' as subset with 1st derivative during the restriction
Vmax<-max(c1) #Identifies the first peak in the 1st derivative curve
TimeIpre<- der1$c.eval.days.[der1$c.nefa.estim.1.==Vmax] #Identifies time when Vmax happens
Ipre<-estimated_data$c.nefa.estim.[estimated_data$c.eval.days.==TimeIpre] #Identifies the inflection point during the restriction

#Define peak during restriction
c2 <- subset(estimated_data, c.eval.days.> 1 & c.eval.days.<= 4, select=c(c.nefa.estim.)) #Creates a subset 'c2' to detect peak
Tt<-max(c2) #Identifies the peak
TimeTt<- estimated_data$c.eval.days.[estimated_data$c.nefa.estim.==Tt] #Identifies time when peak happens

#Determine the end of the recovery period
e2 <- subset(estimated_data, c.eval.days.> 4 & c.eval.days.<= 7, select=c(c.eval.days.,c.nefa.estim.)) #Creates a subset 'e2' with predictions after restriction in order to identify Ppost

#The conditional below calculates TimePpost based on whether or not the smoothed curve crosses baseline
if (e2[which.min(abs(e2$c.nefa.estim.)),]$c.nefa.estim. > Ppre ) {
  TimePpost<-e2[which.min(abs(Ppre-e2$c.nefa.estim.)),]$c.eval.days. #Extracts the moment in time the smoothed curve gets the closest to baseline
  } else {
  TimePpost<-e2[which.min(abs(e2$c.nefa.estim>Ppre)),]$c.eval.days. #Extracts the moment in time the smoothed curve crosses baseline
}

Ppost<-e2[e2$c.eval.days.==TimePpost,]$c.nefa.estim. #Find the value in the smoothed curve that corresponds to TimePpost

#Define inflection point during recovery
c3 <- subset(der1, c.eval.days.> 4 & c.eval.days.<= TimePpost, select=c(c.nefa.estim.1.)) #Creates a subset 'c3' with first derivatives during the recovery, delimited between the end of restriction and TimePpost
Vmin<-min(c3) #Identifies the first peak in the first derivative curve
TimeIpost<- der1$c.eval.days.[der1$c.nefa.estim.1.==Vmin] #Identifies time when Vmin happens
Ipost<-estimated_data$c.nefa.estim.[estimated_data$c.eval.days.==TimeIpost] #Identifies the inflection point during the recovery

#Calculating of timing features
deltaT1<-TimeIpre-TimePpre
deltaT2<-TimeTt-TimePpre
deltaT3<-TimeTt-TimeIpre
deltaT4<-TimeIpost-TimeTt
deltaT5<-TimePpost-TimeTt
deltaT6<-TimePpost-TimeIpost
deltaT7<-TimePpost-TimePpre

Amplitude<-Tt-Ppre
relativeAmplitude=(Amplitude/Ppre)*100
relativeVmin<-Vmin/Amplitude
relativeVmax<-Vmax/Amplitude

#Calculation of areas

#Area under the response curve during the restriction
auc1<-auc(estimated_data$c.eval.days.,estimated_data$c.nefa.estim., from=1, to=4, type = "spline") #Area under curve
aul1<-4*Ppre #Area under baseline
A1<-auc1-aul1

#Area between curves during the recovery period, considering that recovery ends at TimePpost
auc2<-auc(estimated_data$c.eval.days.,estimated_data$c.nefa.estim., from=4, to=TimePpost, type = "spline") #Area under response curve
aul2<-(TimePpost-4)*Ppre #Area under baseline

#Area between curves during the recovery period, considering that recovery ends 3d after the end of restriction
auc3<-auc(estimated_data$c.eval.days.,estimated_data$c.nefa.estim., from=4, to=7, type = "spline") #Area under response curve
aul3<-3*Ppre #Area under baseline

#Conditional statement to calculate A2 according to Timeppost value and theoretical end of challenge
if (TimePpost<=7) {
  A2<-auc2-aul2
  } else {
  A2<-auc3-aul3
}

#Calculation of A3
A3<-A1+A2

#Calculation of recovery variable
recovery<-estimated_data$c.nefa.estim.[estimated_data$c.eval.days.==7] - estimated_data$c.nefa.estim.[estimated_data$c.eval.days.==1]

########################################################################
########## Complete plots ##############################################
########################################################################
nefa<-rep(Ppre, each = 16) 
t<-subset(DF, select=c(t))
baseline<-cbind(t, nefa) #needed to plot baseline


nefa<-subset(DF, select=c(nefa))
t<-rep(0, each = 16) 
beginning_FR<-cbind(t, nefa) #needed to plot beginning of feed restriction


nefa<-subset(DF, select=c(nefa))
t<-rep(4, each = 16) 
end_FR<-cbind(t, nefa) #needed to plot end of feed restriction

nefa<-subset(DF, select=c(nefa))
t<-rep(TimePpost, each = 16) 
end_recovery<-cbind(t, nefa) #needed to plot end of recovery

names(estimated_data)[1]<-paste("t")
names(estimated_data)[2]<-paste("nefa")

plot1 <- ggplot(DF, aes(x=t, y=nefa)) + geom_point(aes(color="Observation")) + geom_line(data = estimated_data, aes(color="Response"))  + 
         geom_line(data = baseline, aes(color="Baseline")) + geom_line(data = beginning_FR, aes(color="Beginning of restriction")) + geom_line(data = end_FR, aes(color="End of restriction")) + geom_line(data = end_recovery, aes(color="End of recovery"))

plot1 +  xlab("Time elapsed (d)") + ylab("Non-esterified fatty acids concentration (g/L)")  + labs(color='Meaning of colours')  + theme_classic()

########################################################################
########## Assemble & export output file 1 #############################
#######################################################################

Variables<-matrix(c("NA","lambda","Ppre","TimePpre","Ipre","TimeIpre","T","TimeT","Ipost","TimeIpost","Ppost","TimePpost","Vmax","Vmin","relativeVmax","relativeVmin","Amplitude","relativeAmplitude","deltaT1","deltaT2","deltaT3","deltaT4","deltaT5","deltaT6","deltaT7","A1","A2","A3","recovery"), nrow=29,ncol=1)
Restriction<-matrix(c("Restriction",lambda,Ppre,TimePpre,Ipre,TimeIpre,Tt,TimeTt,Ipost,TimeIpost,Ppost,TimePpost,Vmax,Vmin,relativeVmax,relativeVmin,Amplitude,relativeAmplitude,deltaT1,deltaT2,deltaT3,deltaT4,deltaT5,deltaT6,deltaT7,A1,A2,A3,recovery),nrow=29,ncol=1)

Output <- cbind(Variables, Restriction)
write.table(Output,file="output.csv")



#Part 2
########################################################################
########## Smoothing MY response #######################################
########################################################################
#This part of the script is adapted to smooth time series and extract dscriptors from the smoothed functions of physiological variables that decrease in face to the perturbation.

cat("\014") #Clear console
rm(list=ls(all=TRUE)) #Remove the history 

DF <- read.table("sample_cow_MY.csv", header = TRUE,sep = ';') 
#This dataset refers to MY response of a cow that was exposed to two feed restrictions. First restriction lasted for 4 days (t=0 d until t = 3 d)
#The dataset has three columns, the first one is time (t) in days (d), the second one is also time (t2, d) and the third column is MY.
#Data file is available for download under https://doi.org/10.5281/zenodo.5776669.

########################################################################
########## Using FDA for MY: baseline ##################################
########################################################################

#Create a sub-set from DF by eliminating variable t
DF2<-subset(DF,select=c(t2,MY))
DF2<- na.omit(DF2)

#Smoothing coefficient and grid parameters:
lambda2 = 1e5 #Penalty coefficient for baseline MY
range2 = c(min(DF2$t2),max(DF2$t2))
break.vec2 = DF2$t2[seq(1, length(DF2$t2),1)]#A vector with the length of the dataset especifying the knots where piecewise polinomials will be fitted.

#Basis function
basis.func2 = create.bspline.basis(rangeval = range2, norder = 6, breaks = break.vec2) #The parameter "norder" is The order of b-splines: cubic, quadratic, etc.
param.func2 = fdPar(fdobj = basis.func2, Lfdobj = 4,lambda = lambda2) #The parameter "Ldfobj" is a linear differential operator object.
smooth.MY2 = smooth.basis(argvals = DF2$t2,DF2$MY,param.func2)
MY.fda2 = smooth.MY2$fd

#Evaluate the estimated function 
eval.days2 = seq(min(DF2$t2),max(DF2$t2),0.1)
MY.estim2   = eval.fd(eval.days2,MY.fda2) #Estimation of MY response at baseline
MY.estim.21 = eval.fd(eval.days2,MY.fda2,1) #Estimation of first derivative from nefa response at baseline

########################################################################
#### Using FDA for MY response during challenges and recoveries ########
########################################################################

#Smoothing coefficient and grid parameters
lambda = 1e3 #Penalty coefficient for response MY
range = c(min(DF$t),max(DF$t))
break.vec = DF$t[seq(1, length(DF$t),1)]#A vector with the length of the dataset especifying the knots where piecewise polinomials will be fitted.

#Basis function
basis.func = create.bspline.basis(rangeval = range, norder = 6, breaks = break.vec) #The parameter "norder" is The order of b-splines: cubic, quadratic, etc.
param.func = fdPar(fdobj = basis.func, Lfdobj = 4,lambda = lambda) #The parameter "Ldfobj" is a linear differential operator object.
smooth.MY = smooth.basis(argvals = DF$t,DF$MY,param.func)
MY.fda = smooth.MY$fd

#Evaluate theestimated function 
eval.days = seq(min(DF$t),max(DF$t),0.1)
MY.estim   = eval.fd(eval.days,MY.fda) #Estimation of MY response
MY.estim.1 = eval.fd(eval.days,MY.fda,1) #Estimation of first derivative from MY response

########################################################################
########## Extracting features from first challenge ####################
########################################################################

#Create tables with information post-FDA for visualisation and further analyses
estimated_data<-data.frame(c(eval.days),c(MY.estim)) 
estimated_baseline<-data.frame(c(eval.days2),c(MY.estim2)) 
der1<-data.frame(c(eval.days),c(MY.estim.1)) 

#Creates a table with estimated data from baseline and response MY
matched<-data.frame(estimated_baseline,estimated_data) #Creates a table with estimated data from baseline and actual MY

#Reference before first challenge
TimePpre1<-0 #Fixed for all cows, set from experimental design
Ppre1<- estimated_baseline$c.MY.estim2.[estimated_baseline$c.eval.days2.==0] 

#Define inflection point during first restriction
#In this example, the first restriction started at t=0 d and ended at t=4 d.
c1 <- subset(der1, c.eval.days.>= 0 & c.eval.days.< 4, select=c(c.MY.estim.1.)) #Creates a subset with first derivatives during the restriction
Vmin1<-min(c1) #Identifies the first valley in the first derivative curve
TimeIpre1<- der1$c.eval.days.[der1$c.MY.estim.1.==Vmin1] #Identifies time when Vmin happens
Ipre1<-estimated_data$c.MY.estim.[estimated_data$c.eval.days.==TimeIpre1] #Identifies the inflection point during the first restriction

#Define peak during first restriction
c2 <- subset(estimated_data, c.eval.days.>= 0 & c.eval.days.< 7, select=c(c.MY.estim.)) #Creates a subset with predictions during the restriction
Tt1<-min(c2) #Identifies the first valley of MY
TimeTt1<- estimated_data$c.eval.days.[estimated_data$c.MY.estim.==Tt1] #Identifies time when first valley happens

#Automatically determine the end of the second recovery period
e2 <- subset(matched, c.eval.days.> 4 & c.eval.days.<= 10, select=c(c.eval.days.,c.MY.estim2., c.MY.estim.)) #Creates a subset with predictions after feed restriction in order to identidy Ppost
e2$diff <- abs(e2$c.MY.estim2. - e2$c.MY.estim. )

#The conditional below calculates TimePpost based on whether or not the smoothed curve crosses baseline
if (e2[which.max(abs(e2$c.MY.estim.)),]$c.MY.estim. < e2[which.min(abs(e2$c.MY.estim2.)),]$c.MY.estim2. ) {
  TimePpost1<-e2[which.max(abs(e2$c.MY.estim.)),]$c.eval.days. #Extracts the moment in time the smoothed curve gets the closest to baseline
  } else {
  TimePpost1<-e2[which.min((e2$diff)),]$c.eval.days. #Extracts the moment in time the smoothed curve crosses baseline
}
Ppost1<-e2[e2$c.eval.days.==TimePpost1,]$c.MY.estim.

#Define inflection point during first recovery
c3 <- subset(der1, c.eval.days.>= 4 & c.eval.days.< TimePpost1, select=c(c.MY.estim.1.)) #Creates a subset with first derivatives during the recovery
Vmax1<-max(c3) #Identifies the first peak in the first derivative curve
TimeIpost1<- der1$c.eval.days.[der1$c.MY.estim.1.==Vmax1] #Identifies time when Vmax happens
Ipost1<-estimated_data$c.MY.estim.[estimated_data$c.eval.days.==TimeIpost1] #Identifies the inflection point during the first recovery

#Calculating of timing features
deltaT11<-TimeIpre1-TimePpre1
deltaT21<-TimeTt1-TimePpre1
deltaT31<-TimeTt1-TimeIpre1
deltaT41<-TimeIpost1-TimeTt1
deltaT51<-TimePpost1-TimeTt1
deltaT61<-TimePpost1-TimeIpost1
deltaT71<-TimePpost1-TimePpre1

#Calculate baseline at Tt1
baselineTt1<-estimated_baseline$c.MY.estim2.[estimated_baseline$c.eval.days2.==TimeTt1]

#Calculation of amplitude from baseline correspondent to the maximum MY
Amplitude1<-baselineTt1-Tt1
relativeVmin1<-Vmin1/Amplitude1
relativeVmax1<-Vmax1/Amplitude1

#Calculating areas

#Area between curves for Challenge 1
auc11<-auc(estimated_data$c.eval.days.,estimated_data$c.MY.estim., from=0, to=4, type = "spline") #Area under red curve
aul11<-auc(estimated_baseline$c.eval.days2.,estimated_baseline$c.MY.estim2., from=0, to=4, type = "spline") #Area under baseline (purple curve)

A11<-aul11-auc11

#Area between curves for recovery 1 considering TimePpost
auc21<-auc(estimated_data$c.eval.days.,estimated_data$c.MY.estim., from=4, to=TimePpost1, type = "spline") #Area under red curve
aul21<-auc(estimated_baseline$c.eval.days2.,estimated_baseline$c.MY.estim2., from=4, to=TimePpost1, type = "spline") #Area under baseline (purple curve)

#Area between curves for recovery 1 considering theoretical ending limit of challenge of 3 days after challenge stopped
auc31<-auc(estimated_data$c.eval.days.,estimated_data$c.MY.estim., from=4, to=7, type = "spline") #Area under red curve
aul31<-auc(estimated_baseline$c.eval.days2.,estimated_baseline$c.MY.estim2., from=4, to=7, type = "spline") #Area under baseline (purple curve)

#Conditional statement to calculate A2 according to Tppost value and theoretical end of challenge
if (TimePpost1<=7) {
  A21<-aul21-auc21
  } else {
  A21<-aul31-auc31
}

#Calculation of A3
A31<-A11+A21

#Calculation of recovery
recovery1<-estimated_data$c.MY.estim.[estimated_data$c.eval.days.==7] - estimated_data$c.MY.estim.[estimated_data$c.eval.days.==0]


########################################################################
########## Extracting features from second challenge ###################
########################################################################

TimePpre2<-21 #Fixed for all cows, set from experimental design
Ppre2<-estimated_baseline$c.MY.estim2.[estimated_baseline$c.eval.days2.==21] #This is determined from "estimated_baseline", when time = 21

#Define inflection point during second restriction
#In this example, the second restriction started at t=21 d and ended at t=31 d.
c4 <- subset(der1, c.eval.days.>= 21 & c.eval.days.< 31, select=c(c.MY.estim.1.)) #Creates a subset with first derivatives during the restriction
Vmin2<-min(c4) #Identifies the second valley in the first derivative curve
TimeIpre2<- der1$c.eval.days.[der1$c.MY.estim.1.==Vmin2] #Identifies time when Vmax happens
Ipre2<-estimated_data$c.MY.estim.[estimated_data$c.eval.days.==TimeIpre2] #Identifies the inflection point during the second restriction

#Define valley during second restriction
c5 <- subset(estimated_data, c.eval.days.>= 21 & c.eval.days.< 31, select=c(c.MY.estim.)) #Creates a subset with predictions during the restriction
Tt2<-min(c5) #Identifies the second valley of MY
TimeTt2<- estimated_data$c.eval.days.[estimated_data$c.MY.estim.==Tt2] #Identifies time when second valley happens

#Automatically determine the end of the second recovery period
e5 <- subset(matched, c.eval.days.> 31 & c.eval.days.<= 40, select=c(c.eval.days.,c.MY.estim2., c.MY.estim.)) #Creates a subset with predictions after feed restriction in order to identidy Ppost
e5$diff <- abs(e5$c.MY.estim2. - e5$c.MY.estim. )

#The conditional below calculates TimePpost based on whether or not the smoothed curve crosses baseline
if (e5[which.max(abs(e5$c.MY.estim.)),]$c.MY.estim. < e5[which.min(abs(e5$c.MY.estim2.)),]$c.MY.estim2. ) {
  TimePpost2<-e5[which.max(abs(e5$c.MY.estim.)),]$c.eval.days. #Extracts the moment in time the smoothed curve gets the closest to baseline
  } else {
  TimePpost2<-e5[which.min((e5$diff)),]$c.eval.days. #Extracts the moment in time the smoothed curve crosses baseline
}

Ppost2<-e5[e5$c.eval.days.==TimePpost2,]$c.MY.estim.

#Define inflection point during second recovery
c6 <- subset(der1, c.eval.days.>= 31 & c.eval.days.< TimePpost2, select=c(c.MY.estim.1.)) #Creates a subset with first derivatives during the recovery
Vmax2<-max(c6) #Identifies the first peak in the first derivative curve
TimeIpost2<- der1$c.eval.days.[der1$c.MY.estim.1.==Vmax2] #Identifies time when Vmax happens
Ipost2<-estimated_data$c.MY.estim.[estimated_data$c.eval.days.==TimeIpost2] #Identifies the inflection point during second recovery

#Calculating of timing features
deltaT12<-TimeIpre2-TimePpre2
deltaT22<-TimeTt2-TimePpre2
deltaT32<-TimeTt2-TimeIpre2
deltaT42<-TimeIpost2-TimeTt2
deltaT52<-TimePpost2-TimeTt2
deltaT62<-TimePpost2-TimeIpost2
deltaT72<-TimePpost2-TimePpre2

#Calculate baseline at Tt2
baselineTt2<-estimated_baseline$c.MY.estim2.[estimated_baseline$c.eval.days2.==TimeTt2]

#Calculation of amplitude from baseline correspondent to the maximum MY
Amplitude2<-baselineTt2-Tt2
relativeVmin2<-Vmin2/Amplitude2
relativeVmax2<-Vmax2/Amplitude2

#Calculating areas

#Area between curves for Challenge 2
auc12<-auc(estimated_data$c.eval.days.,estimated_data$c.MY.estim., from=21, to=31, type = "spline") #Area under red curve
aul12<-auc(estimated_baseline$c.eval.days2.,estimated_baseline$c.MY.estim2., from=21, to=31, type = "spline") #Area under purple curve

A12<-aul12-auc12

#Area between curves for recovery 2
auc22<-auc(estimated_data$c.eval.days.,estimated_data$c.MY.estim., from=31, to=TimePpost2, type = "spline") #Area under red curve
aul22<-auc(estimated_baseline$c.eval.days2.,estimated_baseline$c.MY.estim2., from=31, to=TimePpost2, type = "spline") #Area under purple curve

#Calculation of rebound area after challenge 2
auc32<-auc(estimated_data$c.eval.days.,estimated_data$c.MY.estim., from=31, to=34, type = "spline") #Area under red curve
aul32<-auc(estimated_baseline$c.eval.days2.,estimated_baseline$c.MY.estim2., from=31, to=34, type = "spline") #Area under purple curve

#Conditional statement to calculate A2 according to Tppost value and theoretical end of challenge
if (TimePpost2<=34) {
  A22<-aul22-auc22
  } else {
  A22<-aul32-auc32
}

#Calculation of A3
A32<-A12+A22

#Calculation of recovery
recovery2<-estimated_data$c.MY.estim.[estimated_data$c.eval.days.==34] - estimated_data$c.MY.estim.[estimated_data$c.eval.days.==21]


########################################################################
########## Complete plot ##############################################
########################################################################

MY<-subset(DF, select=c(MY))
t<-rep(0, each = 37) 
beginning_FR1<-cbind(t, MY) #needed to plot beginning of first feed restriction

MY<-subset(DF, select=c(MY))
t<-rep(4, each = 37) 
end_FR1<-cbind(t, MY) #needed to plot end of first feed restriction

MY<-subset(DF, select=c(MY))
t<-rep(TimePpost1, each = 37) 
end_recovery1<-cbind(t, MY) #needed to plot end of first recovery

MY<-subset(DF, select=c(MY))
t<-rep(21, each = 37) 
beginning_FR2<-cbind(t, MY) #needed to plot beginning of second feed restriction

MY<-subset(DF, select=c(MY))
t<-rep(31, each = 37) 
end_FR2<-cbind(t, MY) #needed to plot end of second feed restriction

MY<-subset(DF, select=c(MY))
t<-rep(TimePpost2, each = 37) 
end_recovery2<-cbind(t, MY) #needed to plot end of second recovery

names(estimated_data)[1]<-paste("t")
names(estimated_data)[2]<-paste("MY")

names(estimated_baseline)[1]<-paste("t")
names(estimated_baseline)[2]<-paste("MY")

plot1 <- ggplot(DF, aes(x=t, y=MY)) + geom_point(aes(color="Observation")) + geom_line(data = estimated_data, aes(color="Response")) + geom_line(data = estimated_baseline, aes(color="Baseline")) + 
         geom_line(data = beginning_FR1, aes(color="Beginning of first restriction")) + geom_line(data = end_FR1, aes(color="End of first restriction")) + geom_line(data = end_recovery1, aes(color="End of first recovery")) +
         geom_line(data = beginning_FR2, aes(color="Beginning of second restriction")) + geom_line(data = end_FR2, aes(color="End of second restriction")) + geom_line(data = end_recovery2, aes(color="End of second recovery"))


plot1 +  xlab("Time elapsed (d)") + ylab("Milk yield (kg/d)")  + labs(color='Meaning of colours')  + theme_classic()


########################################################################
########## Assemble & export output file ###############################
########################################################################

Variables<-matrix(c("NA","Ppre","TimePpre","Ipre","TimeIpre","T","TimeT","Ipost","TimeIpost","Ppost","TimePpost","Vmax","Vmin","relativeVmax","relativeVmin","Amplitude","deltaT1","deltaT2","deltaT3","deltaT4","deltaT5","deltaT6","deltaT7","A1","A2","A3","recovery"),nrow=27,ncol=1)
Challenge1<-matrix(c("Challenge 1",Ppre1,TimePpre1,Ipre1,TimeIpre1,Tt1,TimeTt1,Ipost1,TimeIpost1,Ppost1,TimePpost1,Vmax1,Vmin1,relativeVmax1,relativeVmin1,Amplitude1,deltaT11,deltaT21,deltaT31,deltaT41,deltaT51,deltaT61,deltaT71,A11,A21,A31,recovery1),nrow=27,ncol=1)
Challenge2<-matrix(c("Challenge 2",Ppre2,TimePpre2,Ipre2,TimeIpre2,Tt2,TimeTt2,Ipost2,TimeIpost2,Ppost2,TimePpost2,Vmax2,Vmin2,relativeVmax2,relativeVmin2,Amplitude2,deltaT12,deltaT22,deltaT32,deltaT42,deltaT52,deltaT62,deltaT72,A12,A22,A32,recovery2),nrow=27,ncol=1)

Output <- cbind(Variables, Challenge1, Challenge2)
write.table(Output,file="output.csv")