# Prepare R Session -------------------------------------------------------
#sessionInfo()
#R version 4.4.1 (2024-06-14 ucrt)
#Platform: x86_64-w64-mingw32/x64
#Running under: Windows 11 x64 (build 22621)

#set working directory (wd) to where data is stored
setwd ("C:/Users/uk074017/Desktop/KAW_Daten_APs_Ergebnisse_offline/AP1.2 Milchversorgung/Veröffentlichungen/Paper 1/R Project/MSR") ##setwd ("C:/Dokumente/Auswertung")
rm (list= ls ())
if (file.exists ('.RData')) load ('.RData')

#load packages needed for this script
library (triact)
library (blme)
library (DHARMa)
library (boot)
library (pbkrtest)
library(dplyr)
library (vegan)
library (ggplot2)
library (tidyr)
library(forcats)
library(lubridate)

#function for model estimate 
extract.ci <- function (x) {
  out <- data.frame (numeric (0),
                     numeric (0),
                     numeric (0))
  for (i in 1:length (x [['t0']])) {
    out <- rbind (out, c (x [['t0']] [i],
                          boot.ci (x, index= i,
                                   type= 'perc') [['percent']] [, 4:5]))
  }
  names (out) <- c ('estim', 'lo.ci', 'up.ci')
  out
}

# Prepare meta data -------------------------------------------------------

## files and where to find them####################
## list all available files
csv.files <- list.files ('../../MSR-Logger-Data', full.names= TRUE, recursive= TRUE)

## choose files ending with '.csv' of calves
## msr files saved with German names for calf (Kalb)
csv.files <- csv.files [substring (csv.files, nchar (csv.files)-3, nchar (csv.files)) == '.csv']
which.cases <- c (grep ('Kalb', csv.files), grep ('kalb', csv.files))
#which.cases <- which.cases [order (which.cases)]
csv.files <- csv.files [which.cases]

## data.frame with full path and file names
names (csv.files) <- sapply (strsplit (csv.files, '/'), function (x) x [[4]])
names(csv.files) <- basename(csv.files)

##Import Data######################################################################################################
## Before loading here: remove Umlaute from column names in csv file!!
## import other dataframes that might be needed to add missing data via matches
## import MSR data
calves.meta.df <- read.table ('../AP1.3_Calves_MSR.csv',
                              header= TRUE, sep= ';', dec= '.', as.is= FALSE, na.strings= c ('na', '', ' '))

## adjust data types (specifically date-time variable)
calves.meta.df <- cbind (calves.meta.df,
                         data.frame (date.time= strptime (paste (calves.meta.df [, 'Date'],
                                                                 calves.meta.df [, 'Logger_Time']),
                                                          format= '%d.%m.%Y %H:%M:%S', tz= 'GMT')))
calves.meta.df [, 'Phase'] <- factor (as.character (calves.meta.df [, 'Phase']),
                                      levels= c ('1. week', '2. week', '8. week', 'weaning (1)', 'weaning (2)'))
calves.meta.df [, 'CalfID'] <- factor (calves.meta.df [, 'CalfID'])

##Prepare meta/MSR data for triact package##########
## order data by calf ID and time
# omit NN calves that still carry the logger
# calves.meta.df <- calves.meta.df [1:(nrow (calves.meta.df)-NN),]
calves.meta.df <- calves.meta.df [order (calves.meta.df [, 'CalfID'], calves.meta.df [, 'date.time']), ]
summary (calves.meta.df)
dim (calves.meta.df)

## split data set according to logger on / off
calves.on.df <- calves.meta.df [!is.na (calves.meta.df [, 'Logger_new']), ]
calves.off.df <- calves.meta.df [!is.na (calves.meta.df [, 'Logger_old']), ]

## check corrspondence of lines in on / off
dim (calves.on.df)
dim (calves.off.df)

## add NA in case of mismatched numbers of lines
#cbind (rbind (calves.on.df [, c ('date.time', 'CalfID')],
#              data.frame (date.time= rep (NA, 2),
#                          CalfID= rep (NA, 2))),
#       calves.off.df [, c ('date.time', 'CalfID')])
#cbind (calves.on.df [, c ('date.time', 'CalfID')],
#       rbind (calves.off.df [, c ('date.time', 'CalfID')],
#              data.frame (date.time= rep (NA, 2),
#                          CalfID= rep (NA, 2))))

table (as.character (calves.on.df [, 'CalfID']) == as.character (calves.off.df [, 'CalfID']))

## merge needed "on" information to off data
calves.off.df <- cbind (calves.off.df,
                        data.frame (date.time.on= calves.on.df [, 'date.time']))
rm (calves.on.df)

# Create new variables ----------------------------------------------------
## important to have all cells in column date/time formated the same 
calves.off.df [, 'CalWeekMSR'] <- as.numeric (format (calves.off.df [, 'date.time'], '%W'))
calves.off.df [, 'CalWeek']  <- as.numeric (format (calves.off.df [, 'date.time.on'], '%W'))

# checks
table (calves.off.df [, 'CalWeekMSR'] - calves.off.df [, 'CalWeek'])
plot (CalWeekMSR ~ CalWeek, calves.off.df)
abline (a= 1, b= 1)
calves.off.df [calves.off.df [, 'CalWeekMSR'] - calves.off.df [, 'CalWeek'] != 1,]

## add additional information: full path
calves.off.df [, 'full.path'] <- csv.files [as.character(calves.off.df [, 'Filename'])]
calves.off.df [, 'full.path'] <- factor (calves.off.df [, 'full.path'])
rm (csv.files)

# Add additional information: ... ---------------
# e.g. fixed effects: Type of calf = Ammenkalb/eigenes Kalb (Foster calf/ Own calf)
calves.df <- read.table ('../AP1.3_Calves_Overview.csv', fileEncoding = "latin1",
                         header= TRUE, sep= ';', dec= '.', as.is= FALSE, na.strings= c ('na', '', ' '))

calves.df <- calves.df [!is.na(calves.df [,'CalfID']),]
calves.df [duplicated(calves.df [,'CalfID']),]
calves.df[,'CalfID'] <- as.character(calves.df[,'CalfID'])
calves.off.df[,'CalfID'] <- as.character(calves.off.df[,'CalfID'])
row.names (calves.off.df) <- paste (calves.off.df [, 'CalfID'], calves.off.df [, 'CalWeek'])
row.names (calves.df) <- paste (calves.df [, 'CalfID'])

calves.off.df [,'Status'] <-calves.df [calves.off.df[,'CalfID'],'Status']
calves.off.df[, 'LargeGroup'] <- calves.df[calves.off.df[,'CalfID'],'LargeGroup']
calves.off.df [, 'FosterCow'] <- calves.df [calves.off.df[,'CalfID'],'FosterCow'] 
calves.off.df [, 'Breed'] <- calves.df [calves.off.df[,'CalfID'],'Breed'] 
calves.off.df [, 'Sex'] <- calves.df [calves.off.df[,'CalfID'],'Sex'] 

calves.off.df [is.na(calves.off.df [,'Status']),] 

# add location to msr/meta df to match with calender week from calves.health.df
calves.health.df <- read.table ('../AP1.3_Calves_Health.csv', fileEncoding = "latin1",
                             header= TRUE, sep= ';', dec= '.', as.is= FALSE, na.strings= c ('na', '', ' '))
calves.health.df [, 'CalWeek'] <- as.numeric (format (strptime (calves.health.df [, 'Date'],
                                                        format= '%d.%m.%Y', tz= 'GMT'), '%W'))
calves.health.df [, 'location'] <- factor (paste (calves.health.df [, 'Pen_Type'], calves.health.df [, 'Pen_Nr']))

calves.health.df [duplicated (paste (calves.health.df [, 'CalfID'], calves.health.df [, 'CalWeek'])), ]
calves.health.df <- calves.health.df[!is.na(calves.health.df[, 'CalfID']), ]

## Add row names to access information in this table from other tables
row.names(calves.health.df) <- paste (calves.health.df [, 'CalfID'], calves.health.df [, 'CalWeek'])
row.names (calves.off.df) <- paste (calves.off.df [, 'CalfID'], calves.off.df [, 'CalWeek'])

calves.off.df [, 'location'] <- calves.health.df [paste (calves.off.df [, 'CalfID'],
                                                      calves.off.df [, 'CalWeekMSR']), 'location']
calves.off.df [is.na(calves.off.df [,'location']),] 


# ## Check for missing values ---------------------------------------------
## (specifically Phase, file.name - each only once, date.time, date.time.ln
summary (calves.off.df)

calves.off.df [is.na (calves.off.df [, 'Phase']), ]
calves.off.df [is.na (calves.off.df [, 'Filename']), ]
calves.off.df [duplicated (calves.off.df [, 'Filename']), ]    ##NOT FOUND is okay (e.g. lost loggers)
length(calves.off.df [is.na (calves.off.df [, 'full.path']), ])  ##n=23 File not found

str(calves.off.df)
calves.off.df$Date <- dmy(as.character(calves.off.df$Date))

### XXX -------------------------------------------------
## Delete Outliers (1)---------------------------------------------------------
### XXX -------------------------------------------------
#Perfekte Paare
calves_paired <- calves.df %>%
  group_by(FosterCow) %>%
  filter(
    n() == 2,
    all(c("foster calf", "own calf") %in% Status),
    n_distinct(Sex) == 1,
    n_distinct(Breed) == 1
  ) %>%
  mutate(Date_Removal = as.Date(as.character(Date_Removal), format = "%d.%m.%Y")) %>% 
  ungroup()

#Zeitgleiche Daten-> Vergleich zu allen Zeitpunkten
pair_removal <- calves_paired %>%
  group_by(FosterCow) %>%
  summarise(
    Date_Removal_pair = if(all(is.na(Date_Removal))) NA else min(Date_Removal, na.rm = TRUE),
    .groups = "drop"
  )

calves_paired <- calves_paired %>%
  left_join(pair_removal, by = "FosterCow")

#Anwendung auf calves.off.df
calves.off.df <- calves.off.df %>%
  filter(CalfID %in% calves_paired$CalfID)

calves.off.df <- calves.off.df %>%
  mutate(Datum = as.Date(Date)) %>%
  left_join(
    calves_paired %>% select(CalfID, Date_Removal_pair),
    by = "CalfID"
  ) %>%
  filter(
    is.na(Date_Removal_pair) | Datum <= Date_Removal_pair
  ) %>%
  select(-Date_Removal_pair)

nrow(calves.off.df)

#Check: should be 0 (zero)!
calves.off.df %>%
  group_by(FosterCow, Phase) %>%
  summarise(
    n = n(),
    n_FC = sum(Status == "foster calf"),
    n_OC = sum(Status == "own calf"),
    .groups = "drop"
  ) %>%
  filter(n_FC != n_OC )

summary(calves.off.df)
setdiff(unique(calves.off.df$CalfID), unique(calves_paired$CalfID))

# Check: Balance:
calves.off.df %>% 
  group_by(Breed, Sex, Status) %>% 
  summarise(N_CalfIDn = n_distinct(CalfID))

# Check/ Overview plot:
# Structure
summary(calves.off.df$Phase)
str(calves.off.df)
levels(calves.off.df$Phase)
calves.off.df$FosterCow <- as.factor(calves.off.df$FosterCow)
#Plot
ggplot(calves.off.df, aes(x = Phase, y = as.factor(CalfID), color = Status, group = CalfID)) +
  geom_line(linewidth = 0.8, alpha = 0.8) +
  geom_point(size = 2) +
  facet_wrap(~FosterCow, scales = "free_y") +
  labs(
    title = "Sample points of Calves per FosterCow over all Phases",
    subtitle = "Each Facette = One FosterCow, Points show sample points",
    x = "Phase",
    y = "CalfID",
    color = "Status"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    strip.text = element_text(face = "bold"),
    panel.grid.minor = element_blank()
  )


# Which variables to keep -------------------------------------------------
vars.of.interest <- c ('CalfID', 'Batch','FosterCow','LargeGroup', 'location',
                       'date.time.on', 'date.time', 'CalWeek',
                       'Status','Breed','Sex','Phase')

##
### Read accelerometer data (MSR logger) and process with "triact" ------
##

## analysis of MSR-data
options (warn= 1)
rm (msr.calves.df)   #if not run previously, it does not exist (yet)
for (i in (1:nrow(calves.off.df))) {
  #for (i in 1:235) {
  cat ('working on:', 'line', i, '|', as.character (calves.off.df [i, 'full.path']), '\n')
  
  if (!is.na (as.character (calves.off.df [i, 'full.path']))) {
    triact.loc <- Triact$new()
    triact.loc$load_files (input= as.character (calves.off.df [i, 'full.path']),
                           id_substring= c (1, 10),
                           timeFwdUpRight_cols= c(1, -2, 3, -4),
                           skip= '*DATA',
                           dec= '.')
    triact.loc$check_orientation(crit=0.5,interactive=FALSE)        
    triact.loc$add_lying(crit_lie= 0.5, minimum_duration_lying = NULL, window_size = 32)
    #plot(!lying~time,data=triact.loc$data,type="l",bty='n')    
    ints.loc <- triact.loc$summarize_intervals (interval= '24 hours',
                                                duration_units= 'hours',
                                                #lag_in_s= (as.numeric (format (calves.off.df [i, 'date.time.on'],
                                                #                            '%H')) + 2) * 60 * 60,
                                                bouts= TRUE,
                                                calc_for_incomplete= TRUE)
    if (!exists('msr.calves.df')) {
      msr.calves.df <- cbind (calves.off.df [i, vars.of.interest], ints.loc)
    } else {
      msr.calves.df <- rbind (msr.calves.df,
                              cbind (calves.off.df [i, vars.of.interest], ints.loc))
    }
    
    rm (ints.loc)
    rm (triact.loc)
  }
}

rm (i)

msr.calves.df <- msr.calves.df [msr.calves.df [, 'duration'] == 24,]
summary.msr.calves<-summary (msr.calves.df)
dim (msr.calves.df)

table (table (paste (msr.calves.df [, 'CalfID'],  msr.calves.df [, 'CalWeek'])))

### Save after reading MSR ----------
save.image()

### XXX----------
## Delete Outliers (2) -------------------------------------------------------------
### XXX---------------
##If files not found --> pair-partner to be removed too
msr.calves.df %>% 
  group_by(Phase, Status) %>% 
  summarise(N_CalfIDn = n_distinct(CalfID)) %>%  as.data.frame()
table(msr.calves.df$Phase,msr.calves.df$Status)

# Errors:
msr.calves.df %>%
  filter(Phase == "2. week") %>%
  group_by(FosterCow, Status) %>%
  summarise(CalfIDn = list(unique(CalfID)), .groups = "drop") %>%
  pivot_wider(names_from = Status, values_from = CalfIDn) 
#Ammenpaar 84301 Ammenkalb fehlt

# Original rows for pair with inequal rows/data
msr.calves.df %>%
  filter(FosterCow == 84301, Phase == "2. week") %>%
  arrange(date.time.on, CalfID)

msr.calves.df <- msr.calves.df[!(msr.calves.df$Phase == ("2. week") & 
                                   msr.calves.df$FosterCow == 84301), ]

msr.calves.df %>%
  filter(Phase == "weaning(1)") %>%
  group_by(FosterCow, Status) %>%
  summarise(CalfIDn = list(unique(CalfID)), .groups = "drop") %>%
  pivot_wider(names_from = Status, values_from = CalfIDn) 
#In foster cows' pair 32371,50652,84359 foster calf is missing
#In foster cow's pair 758 own calf is missing
msr.calves.df <- msr.calves.df[!(msr.calves.df$Phase == ("weaning(1)") & 
                                   msr.calves.df$FosterCow %in% c(32371,50652,84359,758)), ]

msr.calves.df %>%
  filter(Phase == "weaning (2)") %>%
  group_by(FosterCow, Status) %>%
  summarise(CalfIDn = list(unique(CalfID)), .groups = "drop") %>%
  pivot_wider(names_from = Status, values_from = CalfIDn)
#In foster cow's pair 65801 own calf is missing

msr.calves.df <- msr.calves.df[!(msr.calves.df$Phase == ("weaning (2)") & 
                                   msr.calves.df$FosterCow == 65801), ]

#Check: should be zero rows
msr.calves.df %>%
  group_by(FosterCow, Phase) %>%
  summarise(
    n = n(),
    n_FC = sum(Status == "foster calf"),
    n_OC = sum(Status == "own calf"),
    .groups = "drop"
  ) %>%
  filter(n_FC != n_OC)
# 84306 ok, could be happening due to trimming the msr file (e.g. due to day of removal/loss)
# 32142 and 84359 were equipped according to raw data?

# Original rows for pair with inequal rows/data
msr.calves.df %>%
  filter(FosterCow == 32142, Phase == "1. week") %>%
  arrange(date.time.on, CalfID)

print(
  calves_paired[calves_paired$CalfID %in% msr.calves.df$CalfID, ],
  n = 35
)

msr.calves.df <- msr.calves.df[!(msr.calves.df$Phase == ("1. week") & 
                                   msr.calves.df$FosterCow %in% c(32142, 84359)), ]

#Sample size--------
length(unique(msr.calves.df$CalfID))
n_msr <- 
  msr.calves.df %>% 
  group_by(Phase, Status) %>% 
  summarise(N_CalfIDn = n_distinct(CalfID))
n_msr
table(msr.calves.df$Phase,msr.calves.df$Status)

msr.calves.df %>% 
  group_by(Status) %>% 
  summarise(N_CalfIDn = n_distinct(CalfID))

# Sum contrasts -----------------------------------------------------------
str(msr.calves.df)
contrasts (msr.calves.df [, 'Status']) <- contr.sum (2)
contrasts (msr.calves.df [, 'Phase']) <- contr.sum (5)

msr.calves.df <- cbind(msr.calves.df,
                       model.matrix(~ Status + Phase, data = msr.calves.df) [, -1])
summary (msr.calves.df)
head(msr.calves.df)

## save intermediate step
save.image ()

##
#MODELS----------------------------------------
##

# Model N Bouts Lying -----------------------------------------------------
## model

nBoutsLying.full <- blmer ( log(nBoutsLying) ~
#nBoutsLying.full <- blmer (nBoutsLying ~
                             ## main effects
                             Status1 +
                             Phase1 + Phase2 + Phase3 + Phase4 +
                             ## Two-way interaction
                             Status1:Phase1 + Status1:Phase2 +
                             Status1:Phase3 + Status1:Phase4 +
                             ## Random effect
                            (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                            data = msr.calves.df, 
                            control = lmerControl(optimizer="bobyqa"),
                            REML = FALSE)

summary (nBoutsLying.full)

## model assumptions
nBoutsLying.resid <- simulateResiduals (nBoutsLying.full)
plot (nBoutsLying.resid)
qqnorm (resid (nBoutsLying.full))
scatter.smooth (fitted (nBoutsLying.full ), resid (nBoutsLying.full))
summary(nBoutsLying.resid)

##show outliers
#msr.calves.df[resid (nBoutsLying.full)< -0.5,]
#msr.calves.df[resid (nBoutsLying.full)> 0.5,]
#msr.calves.df[fitted.values (nBoutsLying.full)< 2.0,]
##then check msr curves in msr software (y axis)


plotResiduals (nBoutsLying.resid, msr.calves.df[, 'Status'])
boxplot (split (resid (nBoutsLying.full), msr.calves.df [, 'Status']))

plotResiduals (nBoutsLying.resid, msr.calves.df[, 'Phase'])
boxplot (split (resid (nBoutsLying.full), msr.calves.df [, 'Phase']))

colSums(is.na(msr.calves.df))

qqnorm (ranef (nBoutsLying.full) [['CalWeek']] [, 1])
qqnorm (ranef (nBoutsLying.full) [['location']] [, 1])
qqnorm (ranef (nBoutsLying.full) [['Batch:LargeGroup']] [, 1])
qqnorm (ranef (nBoutsLying.full) [['FosterCow:(Batch:LargeGroup)']] [, 1])
qqnorm (ranef (nBoutsLying.full) [['CalfID:(FosterCow:(Batch:LargeGroup))']] [, 1])
qqnorm (ranef (nBoutsLying.full) [['Phase:(CalfID:(FosterCow:(Batch:LargeGroup)))']] [, 1])

## zero model
nBoutsLying.zero <- blmer (log (nBoutsLying) ~ 1 +
                           (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                           data = msr.calves.df, 
                           control = lmerControl(optimizer="bobyqa"), REML = FALSE)
summary(nBoutsLying.zero)

## single effects
nBoutsLying.status <- blmer (log (nBoutsLying) ~
                           ## main effects
                           #                             Status1 +
                           Phase1 + Phase2 + Phase3 + Phase4 +
                           ## Two-way interaction
                           Status1:Phase1 + Status1:Phase2 +
                           Status1:Phase3 + Status1:Phase4 +
                           ## Random effect
                             (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                           data = msr.calves.df, 
                           control = lmerControl(optimizer="bobyqa"), REML = FALSE)

nBoutsLying.phase <- blmer (log (nBoutsLying) ~
                              ## main effects
                              Status1 +
                              #            Phase1 + Phase1 + Phase3 + Phase4 +
                              ## Two-way interaction
                              Status1:Phase1 + Status1:Phase2 +
                              Status1:Phase3 + Status1:Phase4 +
                              ## Random effect
                              (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                            data = msr.calves.df, 
                            control = lmerControl(optimizer="bobyqa"), REML = FALSE)

nBoutsLying.int <- blmer (log (nBoutsLying) ~
                            ## main effects
                            Status1 +
                            Phase1 + Phase2 + Phase3 + Phase4 +
                            ## Two-way interaction
                            #                              Status:Phase1 + Status:Phase1 +
                            #                              Status:Phase3 + Status:Phase4 +
                            ## Random effect
                            (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                          data = msr.calves.df, 
                          control = lmerControl(optimizer="bobyqa"), REML = FALSE)
##ANOVA-----------
## model comparison (quick but anti-conservative)
# global
anova (nBoutsLying.full, nBoutsLying.zero)

# main effects
anova (nBoutsLying.full, nBoutsLying.status)
anova (nBoutsLying.full, nBoutsLying.phase)

# two-way interaction
anova (nBoutsLying.full, nBoutsLying.int)

## Parametric bootstraps-------
#comparison of models to help decide which is best
# global              
PBmodcomp.nBoutsLying.zero <- PBmodcomp (nBoutsLying.full, nBoutsLying.zero)
summary (PBmodcomp.nBoutsLying.zero)

# main effects
PBmodcomp.nBoutsLying.status <- PBmodcomp (nBoutsLying.full, nBoutsLying.status)
summary (PBmodcomp.nBoutsLying.status)

PBmodcomp.nBoutsLying.phase <- PBmodcomp(nBoutsLying.full, nBoutsLying.phase)
summary (PBmodcomp.nBoutsLying.phase)

## two-way interaction
PBmodcomp.nBoutsLying.int <- PBmodcomp (nBoutsLying.full, nBoutsLying.int)
summary (PBmodcomp.nBoutsLying.int)

## Model estimates ---------------------------------------------------------
estim.df <- data.frame (Status1= rep (c (1, -1), 5),
                        Phase1= rep (c (1, 0, 0, 0, -1), rep(2,5)),
                        Phase2= rep (c (0, 1, 0, 0, -1), rep(2,5)),
                        Phase3= rep (c (0, 0, 1, 0, -1), rep(2,5)),
                        Phase4= rep (c (0, 0, 0, 1, -1), rep(2,5)))

## estimate and confidence interval bootstrap
estim.fn <- function (x) predict (x, newdat= estim.df, re.form= NA)
nBoutsLying.boot <- bootMer (nBoutsLying.full, estim.fn, nsim= 1000, .progress="win")

## back-transformation 
estim.nBoutsLying.df <- cbind (estim.df, exp(extract.ci (nBoutsLying.boot)))

summary(estim.nBoutsLying.df)
estim.nBoutsLying.df


#Estim only status
estim_status.df <- data.frame (Status1= c (1, -1),
                               Phase1=0,
                               Phase2=0,
                               Phase3=0,
                               Phase4=0)

## estimate and confidence interval bootstrap
estim_status.fn <- function (x) predict (x, newdat= estim_status.df, re.form= NA)
nBoutsLying_status.boot <- bootMer (nBoutsLying.full, estim_status.fn, nsim= 1000, .progress="win")

## back-transformation 
estim.nBoutsLying_status.df <- cbind (estim_status.df, exp(extract.ci (nBoutsLying_status.boot)))

summary(estim.nBoutsLying_status.df)
estim.nBoutsLying_status.df

## Graphical presentation ------------------------------------------------------
par(bty = 'n', las = 1) #, mfrow = c (4, 4), mar = c(0.5, 0.5, 0.5, 0.5), oma = c(6, 8, 0, 0))
boxplot (nBoutsLying ~ Status + Phase, msr.calves.df,
         ylab= 'Number lying bouts/ 24h', xlab= '', yaxt= 'n', xaxt= 'n',
         whisklty= 1, range= 0, ylim=c(0,50))
axis (2, seq(0,50,10), seq(0,50,10))
for (i in seq (1, 9, 2)) axis (1, c (i, i+1), c ('FC', 'OC'))
axis (1, seq (1.5, 9.5, 2), c ('week 1', 'week 2', 'week 8', 'pre-\nweaning', 'post-\nweaning'),
      tick= FALSE, line= 2)

for (i in seq (1, 9, 2)) {
  lines (c (i, i+1), estim.nBoutsLying.df [c (i, i+1), 'estim'], lwd= 4)
  lines (c (i, i+1), estim.nBoutsLying.df [c (i, i+1), 'lo.ci'], lwd= 2)
  lines (c (i, i+1), estim.nBoutsLying.df [c (i, i+1), 'up.ci'], lwd= 2)
}
rm (i)


# Model duration lying ----------------------------------------------------
## proportion time lying for logit transformation
## to get values between 0 and 1 (lying % per day)
msr.calves.df [, 'durationLyingLogit'] <- (msr.calves.df [, 'durationLying']/24)

## model
durationLyingLogit.full <- blmer (logit(durationLyingLogit) ~
                              #durationLyingLogit.full <- blmer (durationLyingLogit ~  #as an alternative
                              ## main effects
                              Status1 +
                              Phase1 + Phase2 + Phase3 + Phase4 +
                              ## Two-way interaction
                              Status1:Phase1 + Status1:Phase2 +
                              Status1:Phase3 + Status1:Phase4 +
                              ## Random effect
                                (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                              data = msr.calves.df, 
                              control = lmerControl(optimizer="bobyqa"),
                              REML = FALSE)

summary (durationLyingLogit.full)

## model assumptions
durationLyingLogit.resid <- simulateResiduals (durationLyingLogit.full)
plot (durationLyingLogit.resid)
qqnorm (resid (durationLyingLogit.full))
scatter.smooth (fitted (durationLyingLogit.full ), resid (durationLyingLogit.full))
summary(durationLyingLogit.resid)

##show outliers
#msr.calves.df[resid (durationLyingLogit.full)< -0.6,]
#msr.calves.df[resid (durationLyingLogit.full)> 0.5,]
#msr.calves.df[fitted.values (durationLyingLogit.full)< 0.3,]
#msr.calves.df[fitted.values (durationLyingLogit.full)> 1.8,]
## then check msr curves in msr software (y axis)
## Delete Calf 9126 in 2024-10-01 - 2024-10-08 ? No explanation (sick/lost logger/...) 
#                                                but extreme outlier in fitted.values

plotResiduals (durationLyingLogit.resid, msr.calves.df[, 'Status'])
boxplot (split (resid (durationLyingLogit.full), msr.calves.df [, 'Status']))

plotResiduals (durationLyingLogit.resid, msr.calves.df[, 'Phase'])
boxplot (split (resid (durationLyingLogit.full), msr.calves.df [, 'Phase']))

colSums(is.na(msr.calves.df))

qqnorm (ranef (durationLyingLogit.full) [['CalWeek']] [, 1])
qqnorm (ranef (durationLyingLogit.full) [['location']] [, 1])
qqnorm (ranef (durationLyingLogit.full) [['Batch:LargeGroup']] [, 1])
qqnorm (ranef (durationLyingLogit.full) [['FosterCow:(Batch:LargeGroup)']] [, 1])
qqnorm (ranef (durationLyingLogit.full) [['CalfID:(FosterCow:(Batch:LargeGroup))']] [, 1])
qqnorm (ranef (durationLyingLogit.full) [['Phase:(CalfID:(FosterCow:(Batch:LargeGroup)))']] [, 1])

## zero model
durationLyingLogit.zero <- blmer (logit (durationLyingLogit) ~ 1 +
                             (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                           data = msr.calves.df,
                           control = lmerControl(optimizer="bobyqa"), REML = FALSE)
summary(durationLyingLogit.zero)

## single effects
durationLyingLogit.status <- blmer (logit (durationLyingLogit) ~
                           ## main effects
                           #                             Status1 +
                           Phase1 + Phase2 + Phase3 + Phase4 +
                           ## Two-way interaction
                           Status1:Phase1 + Status1:Phase2 +
                           Status1:Phase3 + Status1:Phase4 +
                           ## Random effect
                             (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                           data = msr.calves.df, 
                           control = lmerControl(optimizer="bobyqa"), REML = FALSE)

durationLyingLogit.phase <- blmer (logit (durationLyingLogit) ~
                              ## main effects
                              Status1 +
                              #                                Phase1 + Phase1 + Phase3 + Phase4 +
                              ## Two-way interaction
                              Status1:Phase1 + Status1:Phase2 +
                              Status1:Phase3 + Status1:Phase4 +
                              ## Random effect
                                (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                              data = msr.calves.df, 
                              control = lmerControl(optimizer="bobyqa"), REML = FALSE)

durationLyingLogit.int <- blmer (logit (durationLyingLogit) ~
                            ## main effects
                            Status1 +
                            Phase1 + Phase2 + Phase3 + Phase4 +
                            ## Two-way interaction
                            #                              Status:Phase1 + Status:Phase1 +
                            #                              Status:Phase3 + Status:Phase4 +
                            ## Random effect
                              (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                            data = msr.calves.df, 
                            control = lmerControl(optimizer="bobyqa"), REML = FALSE)
##ANOVA------------
## model comparison (quick but anti-conservative)
# global
anova (durationLyingLogit.full, durationLyingLogit.zero)

# main effects
anova (durationLyingLogit.full, durationLyingLogit.status)
anova (durationLyingLogit.full, durationLyingLogit.phase)

# two-way interaction
anova (durationLyingLogit.full, durationLyingLogit.int)

## Parametric bootstraps  ----------
#global
PBmodcomp.durationLyingLogit.zero <- PBmodcomp (durationLyingLogit.full, durationLyingLogit.zero)
summary (PBmodcomp.durationLyingLogit.zero)

# main effects
PBmodcomp.durationLyingLogit.status <- PBmodcomp (durationLyingLogit.full, durationLyingLogit.status)
summary (PBmodcomp.durationLyingLogit.status)

PBmodcomp.durationLyingLogit.phase <- PBmodcomp(durationLyingLogit.full, durationLyingLogit.phase)
summary (PBmodcomp.durationLyingLogit.phase)

## two-way interaction
PBmodcomp.durationLyingLogit.int <- PBmodcomp (durationLyingLogit.full, durationLyingLogit.int)
summary (PBmodcomp.durationLyingLogit.int)



## Model estimate ----------------------------------------------------------
estim.df <- data.frame (Status1= rep (c (1, -1), 5),
                        Phase1= rep (c (1, 0, 0, 0, -1), rep(2,5)),
                        Phase2= rep (c (0, 1, 0, 0, -1), rep(2,5)),
                        Phase3= rep (c (0, 0, 1, 0, -1), rep(2,5)),
                        Phase4= rep (c (0, 0, 0, 1, -1), rep(2,5)))

## estimate and confidence interval bootstrap
estim.fn <- function (x) predict (x, newdat= estim.df, re.form= NA)
durationLyingLogit.boot <- bootMer (durationLyingLogit.full, estim.fn, nsim= 1000, .progress="win")

## back-transformation ?
#estim.val <- exp(estim.val)-0.5                 #Counting data
#estim.val <- exp()                              #Amounts
#estim.val <-t(apply(estim.val, 1, inv.logit))   #Proportion
estim.durationLyingLogit.df <- cbind (estim.df, apply (extract.ci (durationLyingLogit.boot), 2, inv.logit))
summary(estim.durationLyingLogit.df)
estim.durationLyingLogit.df


#Estim only status
estim_status.df <- data.frame (Status1= c (1, -1),
                               Phase1=0,
                               Phase2=0,
                               Phase3=0,
                               Phase4=0)

## estimate and confidence interval bootstrap
estim_status.fn <- function (x) predict (x, newdat= estim_status.df, re.form= NA)
durationLyingLogit_status.boot <- bootMer (durationLyingLogit.full, estim_status.fn, nsim= 1000, .progress="win")

## back-transformation ?
#estim.val <- exp(estim.val)-0.5                 #Counting data
#estim.val <- exp()                              #Amounts
#estim.val <-t(apply(estim.val, 1, inv.logit))   #Proportion
estim.durationLyingLogit_status.df <- cbind (estim_status.df, apply (extract.ci (durationLyingLogit_status.boot), 2, inv.logit))
summary(estim.durationLyingLogit_status.df)
estim.durationLyingLogit_status.df

## Graphical presentation---------
par(bty = 'n', las = 1, mar = c(5, 4, 1, 1), oma = c(0.5, 0.5, 0.5, 0.5))
boxplot (durationLyingLogit ~ Status + Phase, msr.calves.df,
         ylab= 'Duration of lying [h/d]', xlab= '', yaxt= 'n', xaxt= 'n',
         whisklty= 1, range= 0, ylim=c(0.25,1))
#axis (2) #, (c(0:1)), c (0:1))
axis (2, seq(0.25,1,0.25), seq(6,24,6))
for (i in seq (1, 9, 2)) axis (1, c (i, i+1), c ('FC', 'OC'))
axis (1, seq (1.5, 9.5, 2), c ('week 1', 'week 2', 'week 8', 'pre-\nweaning', 'post-\nweaning'),
      tick= FALSE, line= 2)
for (i in seq (1, 9, 2)) {
  lines (c (i, i+1), estim.durationLyingLogit.df [c (i, i+1), 'estim'], lwd= 4)
  lines (c (i, i+1), estim.durationLyingLogit.df [c (i, i+1), 'lo.ci'], lwd= 2)
  lines (c (i, i+1), estim.durationLyingLogit.df [c (i, i+1), 'up.ci'], lwd= 2)
}
rm (i)


# Figures GERMAN/DEUTSCH ----------------------------------------------------------
## Graphical presentation
#cols <- rep(c("#4F81BD", "#9BCB59"), times = 5)  # FC blue, OC green
cols <- rep(c("#1F4E79", "#A9D18E"), times = 5)  # FC blue, OC green - more contrast -> Print
col_bw <- rep (c("gray30", "gray80"))

## Figure Lying Bouts--------------
jpeg("../nBoutsLying_PPP.jpg", width = 16, height = 13, units = "cm", res=300, pointsize = 12)
par(bty = 'n', las = 1, mar = c(6, 6, 1, 0.5))
boxplot (nBoutsLying ~ Status + Phase, msr.calves.df,
         ylab= '', xlab= '', yaxt= 'n', xaxt= 'n',
         whisklty= 1, range= 0, ylim=c(0,50),col=cols)
title(ylab = "Anzahl Liegeperioden \npro 24h", cex.lab = 1.2, line=3.0)
axis (2, seq(0,50,10), seq(0,50,10))
for (i in seq (1, 9, 2)) axis (1, c (i, i+1), c ('AK', 'EK'))
axis(1,
     at = seq(1.5, 9.5, 2),
     labels = c("1. LW\n ", "2. LW\n ", "8. LW\n ", "Absetzen\n(1) ", "Absetzen\n(2)"),
     tick = FALSE,
     line = 2.5)
for (i in seq (1, 9, 2)) {
  lines (c (i, i+1), estim.nBoutsLying.df [c (i, i+1), 'estim'], lwd= 4)
  lines (c (i, i+1), estim.nBoutsLying.df [c (i, i+1), 'lo.ci'], lwd= 2)
  lines (c (i, i+1), estim.nBoutsLying.df [c (i, i+1), 'up.ci'], lwd= 2)
}
rm (i)
dev.off()

## Figure Lying Duration--------------
jpeg("../Liegedauer_PPP.jpg", width = 16, height = 13, units = "cm", res=300, pointsize = 12)
par(bty = 'n', las = 1, mar = c(6, 6, 1, 0.5))
boxplot (durationLyingLogit ~ Status + Phase, msr.calves.df,
         ylab= '', xlab= '', yaxt= 'n', xaxt= 'n',
         whisklty= 1, range= 0, ylim=c(0.25,1),col=cols)
#axis (2) #, (c(0:1)), c (0:1))
title(ylab = "Liegedauer \nStunden pro Tag", cex.lab = 1.2, line=3.0)
axis (2, seq(0.25,1,0.25), seq(6,24,6))
for (i in seq (1, 9, 2)) axis (1, c (i, i+1), c ('AK', 'EK'))
axis(1,
     at = seq(1.5, 9.5, 2),
     labels = c("1. LW\n ", "2. LW\n ", "8. LW\n ", "Absetzen\n(1) ", "Absetzen\n(2)"),
     tick = FALSE,
     line = 2.5)
for (i in seq (1, 9, 2)) {
  lines (c (i, i+1), estim.durationLyingLogit.df [c (i, i+1), 'estim'], lwd= 4)
  lines (c (i, i+1), estim.durationLyingLogit.df [c (i, i+1), 'lo.ci'], lwd= 2)
  lines (c (i, i+1), estim.durationLyingLogit.df [c (i, i+1), 'up.ci'], lwd= 2)
}
rm (i)
dev.off()

#Figures ENGLISH ----------------
## Graphical presentation
#cols <- rep(c("#4F81BD", "#9BCB59"), times = 5)  # FC blue, OC green
cols <- rep(c("#1F4E79", "#A9D18E"), times = 5)  # FC blue, OC green - more contrast -> Print
col_bw <- rep (c("gray30", "gray80"))

# Sample size of FC-OC-Pairs per Phase
n.pairs <- tapply(msr.calves.df$FosterCow,
                  msr.calves.df$Phase,
                  function(x) length(unique(na.omit(x))))

## Figure Lying Bouts--------------
jpeg("../nBoutsLying_engl_JDS.jpg", width = 14, height = 10, units = "cm", res=300, pointsize = 8)
par(bty = 'n', las = 1, mar = c(6, 6, 1, 0.5))
bp<- boxplot (nBoutsLying ~ Status + Phase, msr.calves.df,
         ylab= '', xlab= '', yaxt= 'n', xaxt= 'n',
         whisklty= 1, range= 0, ylim=c(0,50),col=cols)
title(ylab = "Number of Lying bouts \nper 24h", cex.lab = 1.2, line=3.0)
axis (2, seq(0,50,10), seq(0,50,10))
for (i in seq (1, 9, 2)) axis (1, c (i, i+1), c ('FC', 'OC'))
axis(1,
     at = seq(1.5, 9.5, 2),
     labels = c("1. week\n ", "2. week\n ", "8. week\n ", "weaning\n(1) ", "weaning\n(2)"),
     tick = FALSE,
     line = 2.5)
# Sample size
axis(1, seq(1.5, 9.5, 2),
     paste0('(n = ', n.pairs, ' pairs)'),
     tick = FALSE, line = 3.5, cex.axis = 0.75)
for (i in seq (1, 9, 2)) {
  lines (c (i, i+1), estim.nBoutsLying.df [c (i, i+1), 'estim'], lwd= 4)
  lines (c (i, i+1), estim.nBoutsLying.df [c (i, i+1), 'lo.ci'], lwd= 2)
  lines (c (i, i+1), estim.nBoutsLying.df [c (i, i+1), 'up.ci'], lwd= 2)
}
rm (i)
dev.off()

## Figure Lying Duration--------------
jpeg("../Liegedauer_englisch_JDS.jpg", width = 14, height = 10, units = "cm", res=300, pointsize = 8)
par(bty = 'n', las = 1, mar = c(6, 6, 1, 0.5))
bp<- boxplot (durationLyingLogit ~ Status + Phase, msr.calves.df,
         ylab= '', xlab= '', yaxt= 'n', xaxt= 'n',
         whisklty= 1, range= 0, ylim=c(0.25,1),col=cols)
#axis (2) #, (c(0:1)), c (0:1))
title(ylab = "Duration of Lying \nper 24h", cex.lab = 1.2, line=3.0)
axis (2, seq(0.25,1,0.25), seq(6,24,6))
for (i in seq (1, 9, 2)) axis (1, c (i, i+1), c ('FC', 'OC'))
axis(1,
     at = seq(1.5, 9.5, 2),
     labels = c("1. week\n ", "2. week\n ", "8. week\n ", "weaning\n(1) ", "weaning\n(2)"),
     tick = FALSE,
     line = 2.5)
# Sample size
axis(1, seq(1.5, 9.5, 2),
     paste0('(n = ', n.pairs, ' pairs)'),
     tick = FALSE, line = 3.5, cex.axis = 0.75)
for (i in seq (1, 9, 2)) {
  lines (c (i, i+1), estim.durationLyingLogit.df [c (i, i+1), 'estim'], lwd= 4)
  lines (c (i, i+1), estim.durationLyingLogit.df [c (i, i+1), 'lo.ci'], lwd= 2)
  lines (c (i, i+1), estim.durationLyingLogit.df [c (i, i+1), 'up.ci'], lwd= 2)
}
rm (i)
dev.off()

# Figures Combined (Same x-Axis) ----------
#cols <- rep(c("#4F81BD", "#9BCB59"), times = 5)  # FC blue, OC green
cols <- rep(c("#1F4E79", "#A9D18E"), times = 5)  # FC blue, OC green - more contrast -> Print
col_bw <- rep (c("gray30", "gray80"))

jpeg("../LyingCombined_engl_dark_JDS.jpg", width = 14, height = 19, units = "cm", res=300, pointsize = 8)
layout(matrix(c(2,1), nrow = 2), heights = c(1,1.5))
par(bty = 'n', las = 1, mar = c(6, 6, 1.5, 0.5), oma=c(0,0,1.5,0))
boxplot (nBoutsLying ~ Status + Phase, msr.calves.df,
              ylab= '', xlab= '', yaxt= 'n', xaxt= 'n',
              whisklty= 1, range= 0, ylim=c(-5,50),col=cols)
title(ylab = "Number of Lying bouts \nper 24h", cex.lab = 1.2, line=3.0)
axis (2, seq(0,50,10), seq(0,50,10))
for (i in seq (1, 9, 2)) axis (1, c (i, i+1), c ('FC', 'OC'))
axis(1,
     at = seq(1.5, 9.5, 2),
     labels = c("1. week\n ", "2. week\n ", "8. week\n ", "weaning\n(1) ", "weaning\n(2)"),
     tick = FALSE,
     line = 2.5)
axis(1, seq(1.5, 9.5, 2),
     paste0('(n = ', n.pairs, ' pairs)'),
     tick = FALSE, line = -2.5, cex.axis = 0.75)
for (i in seq (1, 9, 2)) {
  lines (c (i, i+1), estim.nBoutsLying.df [c (i, i+1), 'estim'], lwd= 4)
  lines (c (i, i+1), estim.nBoutsLying.df [c (i, i+1), 'lo.ci'], lwd= 2)
  lines (c (i, i+1), estim.nBoutsLying.df [c (i, i+1), 'up.ci'], lwd= 2)
}
rm (i)
mtext("b)",
      side = 3, line = 0.5, outer = FALSE,las=0, at=-1, cex=1.5)
# Upper Plot
par(bty = 'n', las = 1, mar = c(1, 6, 1.5, 0.5))

boxplot (durationLyingLogit ~ Status + Phase, msr.calves.df,
              ylab= '', xlab= '', yaxt= 'n', xaxt= 'n',
              whisklty= 1, range= 0, ylim=c(0.25,1),col=cols)
title(ylab = "Duration of Lying \nper 24h", cex.lab = 1.2, line=3.0)
axis (2, seq(0.25,1,0.25), seq(6,24,6))

for (i in seq (1, 9, 2)) {
  lines (c (i, i+1), estim.durationLyingLogit.df [c (i, i+1), 'estim'], lwd= 4)
  lines (c (i, i+1), estim.durationLyingLogit.df [c (i, i+1), 'lo.ci'], lwd= 2)
  lines (c (i, i+1), estim.durationLyingLogit.df [c (i, i+1), 'up.ci'], lwd= 2)
}
rm (i)
mtext("a)",
      side = 3, line = 0.5, outer = FALSE,las=0, at=-1, cex=1.5)
dev.off()

# R verlassen -------------------------------------------------------------
save.image ()
q ('no')

# ___ -------------------------------------------------------------------------