## Prepare R Session -------------------------------------------------------
setwd ("C:/Users/uk074017/Desktop/KAW_Daten_APs_Ergebnisse_offline/AP1.2 Milchversorgung/Veröffentlichungen/Paper 1/R Project/Health") ##setwd ("C:/Dokumente/Auswertung")
rm (list= ls ())
if (file.exists ('.RData')) load ('.RData')



library (triact)
library (blme)
library (DHARMa)
library (boot)
library (pbkrtest)
library (vegan)
library (dplyr)
library (ggplot2)
library (tidyr)
library(forcats)


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
}



##
## Import Data -------------------------------------------------------------
## ProYoungStock Health assessment
health.df <- read.table ('../AP1.3_Calves_Health.csv',
                        header= TRUE, sep= ';', dec= ',', as.is= FALSE, fill = TRUE, na.strings= c ('na', '', ' '))
health.df [, 'CalfID'] <- as.character (health.df [, 'CalfID'])
health.df$Weight <- as.character(health.df$Weight)
health.df$Weight <- gsub(",", ".", health.df$Weight)
summary(health.df$Weight)
health.df$Weight <- as.numeric(health.df$Weight) #. as decimal! in csv file
health.df [, 'CalWeek'] <- as.numeric (format (strptime (health.df [, 'Date'],
                                                    format= '%d.%m.%Y', tz= 'GMT'), '%W'))
health.df [, 'Date'] <- as.Date(health.df[, 'Date'], format= '%d.%m.%Y')
health.df [, 'location'] <- factor (paste (health.df [, 'Pen_Type'], health.df [, 'Pen_Nr']))
health.df [, 'Phase'] <- factor (as.character (health.df [, 'Phase']),
                                    levels= c ('small group', 'large group (beginning)',
                                               'large group (mid)', 'large group (end)',
                                               'weaning (1)', 'weaning (2)',
                                               'weaned'))
summary (health.df)
health.df<-health.df [!is.na(health.df [,'CalfID']),]

# Include further variables (fixed and random effects)----------
calves.df <- read.table ('../AP1.3_Calves_Overview.csv', fileEncoding = "latin1",
                         header= TRUE, sep= ';', dec= '.', as.is= FALSE, na.strings= c ('na', '', ' '))
calves.df [, 'CalfID'] <- as.character (calves.df [, 'CalfID'])
calves.df [, 'Birthdate'] <- as.Date(calves.df[, 'Birthdate'], format= '%d.%m.%Y')
calves.df <- calves.df [!is.na(calves.df [,'CalfID']),]
calves.df [duplicated(calves.df [,'CalfID']),]

row.names (calves.df) <- as.character (calves.df [, 'CalfID'])

health.df [, 'Status'] <- calves.df [health.df[,'CalfID'],'Status']  
health.df [, 'Birthdate'] <- calves.df [health.df[,'CalfID'],'Birthdate']
health.df[, 'LargeGroup'] <- calves.df[health.df[,'CalfID'],'LargeGroup']
health.df [, 'FosterCow'] <- as.factor(calves.df [health.df[,'CalfID'],'FosterCow'])  
health.df [, 'Breed'] <- as.factor(calves.df [health.df[,'CalfID'],'Breed'])  
health.df [, 'Sex'] <- as.factor(calves.df [health.df[,'CalfID'],'Sex'])  

health.df[, 'Age'] <- floor(as.numeric(difftime(health.df[, 'Date'], health.df[, 'Birthdate'], units = "weeks")))
health.df[, 'Age_Days'] <- floor(as.numeric(difftime(health.df[, 'Date'], health.df[, 'Birthdate'], units = "days")))

# Phasen Zuordnung automatisiert--> Phase_Test -------------------------------
# Sicherstellen, dass Date ein Date-Objekt ist
health.df$Date <- as.Date(health.df$Date)

# Part 1: Assigning based on location
health.df$Phase_Test <- NA  # Neue Spalte anlegen
# Small group
health.df$Phase_Test[grepl("small group", health.df$location)] <- "small group"
# Weaning
health.df$Phase_Test[health.df$location == "weaning pen 6"] <- "weaning (1)"
health.df$Phase_Test[health.df$location == "weaning pen 4"] <- "weaning (2)"
health.df$Phase_Test[health.df$location %in% c("weaning pen 1", "weaning pen 2", "weaning pen 3")] <- "weaned"
# Intermediate step: large group
largegroup_locations <- c("large group 2", "large group 3",
                            "large group 4", "large group 5", "large group 6",
                            "large group 7", "large group 8", "large group 9")
health.df$Phase_Test[health.df$location %in% largegroup_locations] <- "large group"

# Part 2: Timely subsets in large groups 
# subset for large group
largegroup.df <- health.df[health.df$Phase_Test == "large group", ]
# Start into large group based on FosterCow
duration <- largegroup.df %>%
  group_by(FosterCow) %>%
  summarise(
    start = min(Date),
    end = max(Date),
    .groups = "drop"
  )

largegroup.df <- largegroup.df %>%
  left_join(duration, by = "FosterCow") %>%
  mutate(
    days_since_start = as.numeric(Date - start),
    Phase_Test = case_when(
      days_since_start <= 21 | location %in% c("large group 2", "large group 3") ~ "large group (beginning)",
      location %in% c("large group 7", "large group 8", "large group 9")         ~ "large group (end)",
      TRUE                                                                      ~ "large group (mid)"
    )
  )

# Part 3: Replace old large group rows in health.df
health.wO.largegroup <- health.df[health.df$Phase_Test != "large group", ]
# Kombinieren
largegroup.df <- largegroup.df[, names(health.df)]
health.df <- rbind(health.wO.largegroup, largegroup.df)
#Vergleich/Kontrolle
table(health.df$Phase, health.df$Phase_Test)
table(health.df$Phase_Test, health.df$location)
table(health.df$Phase, health.df$location)
#Phase_Test nutzen:
health.df$Phase <- NULL
names(health.df)[names(health.df) == "Phase_Test"] <- "Phase"
health.df$Phase <- as.factor(health.df$Phase)

health.df [, 'Phase'] <- factor (as.character (health.df [, 'Phase']),
                                 levels= c ('small group', 'large group (beginning)',
                                            'large group (mid)', 'large group (end)',
                                            'weaning (1)', 'weaning (2)',
                                            'weaned'))
summary(health.df)
#seltsamerweise sind 3 NA Zeilen da, deswegen:
health.df<-health.df [!is.na(health.df [,'CalfID']),] 

## Zuordnungs Check: ob alles aus calves.df richtig zugeordnet wurde-----------
names(health.df)
health.check.df <- (health.df [, c (18,1,21,16,22,6,2,3,19,5,20,23,27)])
summary (health.check.df)
health.df [is.na(health.df [,'Phase']),]
health.df [is.na(health.df [,'Status']),]
health.df [is.na(health.df [,'location']),] 
health.df [is.na(health.df [,'Birthdate']),]
health.df [is.na(health.df [,'FosterCow']),]
health.df [is.na(health.df [,'LargeGroup']),]

##Alle zufälligen Effekte müssen as.factors sein!
health.df[,'CalfID'] <- as.factor(health.df[,'CalfID'])
health.df[,'LargeGroup'] <- as.factor(health.df[,'LargeGroup'])
health.df[,'FosterCow'] <- as.factor(health.df[,'FosterCow'])
health.df[,'location'] <- as.factor(health.df[,'location'])
health.df[,'CalWeek'] <- as.factor(health.df[,'CalWeek'])

# Score -------------------------------------------------------------
health.df[, 'Score'] <- (health.df [, 'Vitality']+
                                 health.df [, 'BCS']+
                                 health.df [, 'Coat_Condition']+
                                 health.df [, 'Diarrhea']+
                                 health.df [, 'Coughing']+
                                 health.df [, 'Ocular_Discharge']+
                                 health.df [, 'Nasal_Discharge']+
                                 health.df [, 'Breathing'])
# These are the exact PYS (Ivemeyer et al. 2022) indikators
# Also,'Hairloss' was not recorded from the beginning and results in NAs in Score

health.df$Score <- as.numeric(health.df$Score)
health.df <- subset(health.df, !(is.na(Weight) & is.na(Score)))

summary(health.df$Score)
#32 NAs vorwiegend beim Absetzen

### XXX-------------------------------------------------
## Delete Outliers (1) ---------------------------------------------------------
### XXX -------------------------------------------------

#Perfect pairs
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()

#Equal Data-> Comparison of BOTH at all times
removal_pair <- 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(removal_pair, by = "FosterCow")

#Applyin to health.df
health.df <- health.df %>%
  filter(CalfID %in% calves_paired$CalfID)

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

nrow(health.df)
#Should be an even number
#Needs removal of individual rows?

#health.df <- health.df %>%
#  group_by(FosterCow, Date) %>%
#  # markiere Partner, deren Gegenstück am selben Date im Krankenstall ist
#  mutate(partner_im_KS = any(Boxtyp == "Krankenstall") & Boxtyp != "Krankenstall") %>%
#  # lösche nur diese Partner-Zeilen
#  filter(!(partner_im_KS & Boxtyp != "Krankenstall")) %>%
#  select(-partner_im_KS) %>%
#  ungroup()

#Cases in which the partner of a pair is missing
tobe_removed <- data.frame(
  Date = as.Date(c("2024-06-24", "2024-06-11", "2024-07-16")),
  CalfID = c(38703, 9147, 38748))

health.df <- health.df[ !paste(health.df$Date, health.df$CalfID) %in% 
                          paste(tobe_removed$Date, tobe_removed$CalfID), ]


#Kontrolle (sollte 0 Zeilen sein)
health.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 )



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

calves_not_paired <- calves.df[ 
  !calves.df$CalfID %in% calves_paired$CalfID, 
]
summary(calves_not_paired)
save.image()

# Calculate Average Daily Weight gain-------------------------------------
names(health.df)
weight.df <- health.df[,c(1:6,18:28)]
weight.df <- weight.df[order(weight.df$CalfID, weight.df$Age_Days), ]

weight.df[is.na(weight.df$Weight),]         #Only in last phases 
length(weight.df[is.na(weight.df$Weight),]) #n=17
weight.df <- weight.df[!is.na(weight.df$Weight),]
summary(weight.df)                          #never weight when fully weaned (only health assessment here) 
table(weight.df$Phase, weight.df$Status)

weight.df$Weight_diff <- c(NA, diff(weight.df$Weight))
weight.df$Days_diff <- c(NA, diff(weight.df$Date))

weight.df$rADWG <- NA
for (i in 2:nrow(weight.df)) {
  weight.df$rADWG[i] <- 
    (weight.df$Weight[i] - weight.df$Weight[i-1]) / weight.df$Weight[i-1]/ weight.df$Days_diff[i]
}
weight.df$ADWG <- weight.df$Weight_diff / weight.df$Days_diff

weight.df[is.na(weight.df$rADWG),] #weil für erste Wiegung gibts keinen Bezug
weight.df[is.na(weight.df$ADWG),] #weil für erste Wiegung gibts keinen Bezug
weight.df <- weight.df[!(weight.df$FosterCow == 84252 & weight.df$Date == as.Date("2024-04-29")), ]


weight.df[weight.df$Days_diff<0,]             #reference is the row (other calf) above
weight.df<-weight.df[weight.df$Days_diff>0,]
summary(weight.df)
weight.df <- weight.df[!is.na(weight.df),]
weight.df <- weight.df[!is.na(weight.df$Phase),]
summary(weight.df)

###XXX-------------------------------------------------
## Delete Outliers (2) -----------------------------
###XXX-------------------------------------------------

### Cleaning from ADWG-Calculation----------------
#Removal of unrealistic ADWG
hist(weight.df$ADWG,breaks=50)
weight.df[weight.df$ADWG<(-1),]
#Outlier+Partner removed:
weight.df <- weight.df[!(weight.df$FosterCow == 50520 & weight.df$Date == as.Date("2024-06-24")), ]

weight.df[weight.df$ADWG>2,]

#Number of calves with weight loss:
length(weight.df[weight.df$ADWG<0,]) #21
summary(weight.df[weight.df$ADWG<0,])
weight.df[weight.df$ADWG<0,]

boxplot (ADWG ~ Status + Phase, weight.df)

#Show outliers
outlier_rows <- do.call(
  rbind,
  by(weight.df, list(weight.df$Status, weight.df$Phase), function(d) {
    o <- boxplot.stats(d$ADWG)$out
    d[d$ADWG %in% o, ]
  })
)

outlier_rows


length(unique(weight.df$CalfID))
n_weight <-
  weight.df %>% 
  group_by(Phase, Status) %>% 
  summarise(N_CalfIDn = n_distinct(CalfID))
n_weight
table(weight.df$Phase,weight.df$Status)


# Check, how CalfIDn per Phase are unbalanced
weight.df %>%
  filter(Phase=="small group") %>% 
  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 )

# Show original rows for pair:
weight.df %>%
  filter(FosterCow == 32566, Phase == "small group") %>%
  arrange(Date, CalfID)
#why not a problem in health.df?
health.df %>%
  filter(FosterCow == 32566, Phase == "small group") %>%
  arrange(Date, CalfID)
#Weighing mistake with partner on this day
weight.df <- weight.df[!(weight.df$Date == as.Date("2024-05-28") & 
                           weight.df$CalfID == 9133), ]


# Reihenfolge der Phasen beibehalten
phases_order <- c ('small group', 'large group (beginning)',
                   'large group (mid)', 'large group (end)',
                   'weaning (1)', 'weaning (2)',
                   'weaned')

weight.df <- weight.df %>%
  mutate(
    Phase = factor(Phase, levels = phases_order),
    FosterCow = as.factor(FosterCow))

ggplot(weight.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 = "Timing per calf pair/foster cow per phase",
    subtitle = "Each Facette = 1 FosterCow, Points show phases",
    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()
  )

table(weight.df$Days_diff)
weight.df[weight.df$Days_diff>35,]


# LifeNetGain-------------
LifeNetGain <- weight.df %>%
  group_by(CalfID) %>%
  arrange(Date) %>%
  mutate(
    NetGain = (last(Weight) - first(Weight)) / last(Age_Days)
  ) %>%
  slice_tail(n = 1) %>%
  ungroup()

boxplot(NetGain ~ Status, LifeNetGain)

summary(LifeNetGain)

LifeNetGain %>% 
  group_by(Status) %>% 
  summarise(NetGain=mean(NetGain))

weight.df %>% 
  group_by(Status) %>% 
  summarise(NetGain=mean(ADWG))
weight.df %>% 
  group_by(Phase, Status) %>% 
  summarise(NetGain=mean(ADWG))

# Final Weight  ------------------------------------------------------
final.weight.df <- health.df %>%
  filter(Phase=='weaning (2)') %>% 
  group_by(CalfID) %>%
  arrange(desc(Date)) %>%
  slice(1)                              

table(final.weight.df$Phase,final.weight.df$Age)

names(final.weight.df)
final.weight.df <- final.weight.df[,c(1:6,16,18:28)]
summary(final.weight.df)

final.weight.df <- droplevels(final.weight.df[final.weight.df$Phase == 'weaning (2)', ])

summary(final.weight.df)
final.weight.df<-final.weight.df[!is.na(final.weight.df$Weight),]
summary(final.weight.df)
final.weight.df <- as.data.frame(final.weight.df)

#Impact/Influence of Age
table(final.weight.df$Age)
table(final.weight.df$Age_Days)
hist(final.weight.df$Age, breaks=5)
hist(final.weight.df$Age_Days, breaks=10)
# According to literature, this has little influence
# No need to include in the model

# Some animals left the sample earlier
# sample reduced until "weaning (2)"
length(unique(final.weight.df$CalfID)) #n=28 Tiere für Nettozunahmen
summary(final.weight.df)

## Intermediate Save ------------------
save.image()

# # Prävalenzberechnung -------------------------------------------------------
summary_health_status <- aggregate(cbind(Vitality,
                                         BCS,Coat_Condition,Diarrhea, Coughing, Nasal_Discharge,
                                         Ocular_Discharge,Breathing) 
                                   ~ Status, data = health.df, FUN = function(x) round(mean(x == 1, na.rm = TRUE)*100, digits=2))
print(summary_health_status) 

summary_health_status_phase <- aggregate(cbind(Vitality,
                                               BCS,Coat_Condition,Diarrhea, Coughing, Nasal_Discharge,
                                               Ocular_Discharge,Breathing) 
                                         ~ Status + Phase, 
                                         data = health.df, 
                                         FUN = function(x) round(mean(x == 1, na.rm = TRUE)*100, digits = 2))

# n per Phase-Status-Group berechnen
n_per_group <- aggregate(Vitality ~ Status + Phase, 
                         data = health.df, 
                         FUN = function(x) sum(!is.na(x)))
colnames(n_per_group)[3] <- "n"

summary_health_status_phase <- merge(summary_health_status_phase, n_per_group, 
                                     by = c("Phase", "Status"))
summary_health_status_phase <- summary_health_status_phase[
  order(factor(summary_health_status_phase$Phase,
               levels = levels(health.df$Phase))),]


print(summary_health_status_phase)
write.csv(summary_health_status_phase, file = "Prevalence_Status_and_Phase.csv", row.names = FALSE)

# PCA ---------------------------------------------------------------------
#decided against PCA, and pro Score as model

#Preparation: pca cannot be calculated if NA are present
#health.df <- na.omit(health.df)
#summary(health.df)

# PCA inklusive PYS, Weight und rADWG
#health.pca <- princomp (health.df[,c (5,22,6,7:14)], cor= TRUE)
#summary (health.pca)
#eigenvals (health.pca)
#plot(health.pca)
# print (loadings (health.pca), cutoff= 0)
#print (loadings (health.pca), cutoff= 0.3)
# print (loadings (health.pca), cutoff= 0.4)
# print (loadings (health.pca), cutoff= 0.5)

#health.comp.df <- rbind (predict (health.pca) [, 1:5]) # for the first X PCs
#health.df <- cbind(health.df, health.comp.df)
#summary(health.df)



# Sample size -------------------------------------------------------
# Check: Balance:
Balance.df <- health.df %>% 
  group_by(Breed, Sex, Status) %>% 
  summarise(N_CalfIDn = n_distinct(CalfID))

ggplot(Balance.df, aes(x = Sex, y = N_CalfIDn, fill = Status)) +
  geom_col(position = position_dodge(width = 0.8)) +
  facet_wrap(~ Breed) +   # Trennt nach Sex
  labs(
    x = "Sex",
    y = "Anzahl CalfIDn",
    fill = "Status"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  )


ggplot(Balance.df, aes(x = Status, y = N_CalfIDn, fill = Sex)) +
  geom_col(position = position_dodge(width = 0.8)) +
  facet_wrap(~ Breed) +   # Trennt nach Sex
  labs(
    x = "Status",
    y = "Anzahl CalfIDn",
    fill = "Status"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

length(unique(health.df$CalfID))
n_health <- 
  health.df %>% 
  group_by(Phase, Status) %>% 
  summarise(N_CalfIDn = n_distinct(CalfID))
n_health
table(health.df$Phase,health.df$Status)

ggplot(n_health, aes(x = Phase, y = N_CalfIDn, group = Status, color = Status)) +
  geom_col(aes(fill = Status), position = "dodge", alpha = 0.6) +
  geom_line(position = position_dodge(width = 0.9), linewidth = 1) +
  geom_point(position = position_dodge(width = 0.9), size = 2) +
  labs(title = "Number of Animals per Phase and Status",
       x = "Phase", y = "Number of Animals") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

length(unique(weight.df$CalfID))
weight.df %>% 
  group_by(Breed, Sex, Status) %>% 
  summarise(N_CalfIDn = n_distinct(CalfID))

n_weight <-
  weight.df %>% 
  group_by(Phase, Status) %>% 
  summarise(N_CalfIDn = n_distinct(CalfID))
n_weight
table(weight.df$Phase,weight.df$Status)
ggplot(n_weight, aes(x = Phase, y = N_CalfIDn, group = Status, color = Status)) +
  geom_col(aes(fill = Status), position = "dodge", alpha = 0.6) +
  geom_line(position = position_dodge(width = 0.9), linewidth = 1) +
  geom_point(position = position_dodge(width = 0.9), size = 2) +
  labs(title = "Weight: Anzahl Tiere pro Phase nach Status",
       x = "Phase", y = "Anzahl Tiere") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
# When calves were in the small group only for 1 week, then no value, 
# Thats why sample size increases for large group beginning

#Final sample
final.weight.df %>% 
  group_by(Status, Breed, Sex) %>% 
  summarise(N_CalfIDn = n_distinct(CalfID))

final.weight.df %>% 
    group_by(Breed, Sex, Status) %>% 
    summarise(Mean_Weight= mean(Weight),
              Median_Weight=median(Weight),
              Min_Weight= min(Weight),
              Max_Weight= max(Weight),
              SD_Weight= sd(Weight)
              )

length(unique(final.weight.df$CalfID))

## after dplyr formating from tibble into dataframe needed
health.df <- as.data.frame(health.df)
weight.df <- as.data.frame(weight.df)
final.weight.df <- as.data.frame(final.weight.df)
calves.df <- as.data.frame(calves.df)


# Models -----------------------------------------------------------------
## Sum contrasts -----------------------------------------------------------
# health.df
str(health.df)
contrasts (health.df [, 'Status']) <- contr.sum (2)
length(unique(health.df$Phase))
contrasts (health.df [, 'Phase']) <- contr.sum (7)
contrasts (health.df [, 'Sex']) <- contr.sum (2)

health.df <- cbind (health.df,
                    model.matrix(~ Status + Phase+ Sex, data = health.df) [, -1])

summary (health.df)
head(health.df)

#weight.df
str(weight.df)
length(unique(weight.df$Phase))
levels(weight.df$Phase)
weight.df <- weight.df[weight.df$Phase!="weaned",]
weight.df <- droplevels(weight.df[weight.df$Phase != "weaned", ])

contrasts (weight.df [, 'Status']) <- contr.sum (2)
contrasts (weight.df [, 'Phase']) <- contr.sum (6)
contrasts (weight.df [, 'Sex']) <- contr.sum (2)

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

#final.weight.df
contrasts (final.weight.df [, 'Status']) <- contr.sum (2)
contrasts (final.weight.df [, 'Sex']) <- contr.sum (2)

final.weight.df <- cbind (final.weight.df,
                    model.matrix(~ Status+ Sex, data = final.weight.df) [, -1])

#contrasts (final.weight.df [, 'Age']) <-  #as continous Covariate/Variabel tranformed (tf)
final.weight.df$Age_tf <- scale(final.weight.df$Age, center = TRUE, scale = TRUE)
summary(final.weight.df$Age_tf)         #Age_tf stands for transformed
summary (final.weight.df)
head(final.weight.df)
names(final.weight.df)

## save intermediate step
save.image()


# Model Score Gesundheit --------------------------------------------------
### Log Transformation ----------------------------------------------------
# log(x+0.5) transformation weil ganzzahlige scores, die sich vielleicht so ähnlich wie Zähldaten verhalten
# Log-Transformation nur auf die numerischen Spalten
# health.df[ , sapply(health.df, is.numeric)] <- log(health.df[ , sapply(health.df, is.numeric)] + 0.5)
# Nicht möglich/nötig bei  0/1 Daten!!

health.df [is.na(health.df [,'Score']),]
health.df <- health.df [!is.na(health.df [,'Score']),]    

## log transformation
health.df$GesamtscoreLog <- log(health.df$Score +0.5)

## Bestes Modell und Optimizer! Poisson, bglmer, binomial etc alles ausprobiert!
##Log model
GesamtscoreLog.full <- blmer (GesamtscoreLog ~
                                ## main effects
                                Status1 +
                                Phase1 + Phase2 + Phase3 + Phase4 +
                                Phase5 + Phase6 +
                                ## Two-way interaction
                                Status1:Phase1 + Status1:Phase2 +
                                Status1:Phase3 + Status1:Phase4 +
                                Status1:Phase5 + Status1:Phase6 +
                                ## Random effect
                                (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                              data = health.df, REML = FALSE)

summary (GesamtscoreLog.full)

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

#Show outliers
#health.df[resid (GesamtscoreLog.full)>  1.5,]
#health.df[resid (GesamtscoreLog.full)< -0.6,]
#health.df[fitted.values (GesamtscoreLog.full)< -2.0,]

plotResiduals (GesamtscoreLog.resid, health.df[, 'Status'])
boxplot (split (resid (GesamtscoreLog.full), health.df [, 'Status']))

plotResiduals (GesamtscoreLog.resid, health.df[, 'Phase'])            
boxplot (split (resid (GesamtscoreLog.full), health.df [, 'Phase']))

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

ranef (GesamtscoreLog.full) [['FosterCow:(Batch:LargeGroup)']] 
##Outlier: 32368:B:1
qqnorm (ranef (GesamtscoreLog.full) [['CalfID:(FosterCow:(Batch:LargeGroup))']] [, 1])
ranef (GesamtscoreLog.full) [['CalfID:(FosterCow:(Batch:LargeGroup))']] 
#Outlier: 38720:32334:E:2
qqnorm (ranef (GesamtscoreLog.full) [['Phase:(CalfID:(FosterCow:(Batch:LargeGroup)))']] [, 1])
ranef (GesamtscoreLog.full) [['Phase:(CalfID:(FosterCow:(Batch:LargeGroup)))']] 
#Outliers:
# GrossBatch_Anfang:38749:50652:G:3
# GrossBatch_Anfang:38720:32334:E:2
# GrossBatch_Anfang:38687:32554:B:1


## zero model
GesamtscoreLog.zero <- blmer (GesamtscoreLog ~ 1 +
                                ## Random effect
                                (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                              data = health.df, REML = FALSE)
summary(GesamtscoreLog.zero)

## single effects
GesamtscoreLog.status <- blmer (GesamtscoreLog ~
                                  ## main effects
                                  #                        Status1 +
                                  Phase1 + Phase2 + Phase3 + Phase4 +
                                  Phase5 + Phase6 +
                                  ## Two-way interaction
                                  Status1:Phase1 + Status1:Phase2 +
                                  Status1:Phase3 + Status1:Phase4 +
                                  Status1:Phase5 + Status1:Phase6 +
                                  ## Random effect
                                  (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                                data = health.df, REML = FALSE)

GesamtscoreLog.phase <- blmer (GesamtscoreLog ~
                                 ## main effects
                                 Status1 +
                                 #                           Phase1 + Phase2 + Phase3 + Phase4 + ...
                                 ## Two-way interaction
                                 Status1:Phase1 + Status1:Phase2 +
                                 Status1:Phase3 + Status1:Phase4 +
                                 Status1:Phase5 + Status1:Phase6 +
                                 ## Random effect
                                 (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                               data = health.df,REML = FALSE)

GesamtscoreLog.int <- blmer (GesamtscoreLog ~
                               ## main effects
                               Status1 +
                               Phase1 + Phase2 + Phase3 + Phase4 +
                               Phase5 + Phase6 +
                               ## Two-way interaction
                               #                         Status1:Phase1 + Status1:Phase2 +
                               #                         Status1:Phase3 + Status1:Phase4 +
                               ## Random effect
                               (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                             data = health.df, REML = FALSE)

### ANOVA model comparison (quick but anti-conservative)---------------
# global
anova (GesamtscoreLog.full, GesamtscoreLog.zero)

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

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

save.image()

### Bootstrap GesamtscoreLog --------------------------------------------------------
# global
PBmodcomp.GesamtscoreLog.zero <- PBmodcomp (GesamtscoreLog.full, GesamtscoreLog.zero)

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

PBmodcomp.GesamtscoreLog.phase <- PBmodcomp(GesamtscoreLog.full, GesamtscoreLog.phase)

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

#Summaries
summary (PBmodcomp.GesamtscoreLog.zero)
summary (PBmodcomp.GesamtscoreLog.status)
summary (PBmodcomp.GesamtscoreLog.phase)
summary (PBmodcomp.GesamtscoreLog.int)


### Model estimate ------------
estim.df <- data.frame (Status1= rep (c (1, -1), 7),
                        Phase1= rep (c (1, 0, 0, 0, 0, 0, -1), rep(2,7)),
                        Phase2= rep (c (0, 1, 0, 0, 0, 0, -1), rep(2,7)),
                        Phase3= rep (c (0, 0, 1, 0, 0, 0, -1), rep(2,7)),
                        Phase4= rep (c (0, 0, 0, 1, 0, 0, -1), rep(2,7)),
                        Phase5= rep (c (0, 0, 0, 0, 1, 0, -1), rep(2,7)),
                        Phase6= rep (c (0, 0, 0, 0, 0, 1, -1), rep(2,7)),
                        Phase7= rep (c (0, 0, 0, 0, 0, 0, -1), rep(2,7)),
                        Sex1 = 0)
## estimate and confidence interval bootstrap
estim.fn <- function (x) predict (x, newdat= estim.df, re.form= NA)
GesamtscoreLog.boot <- bootMer (GesamtscoreLog.full, estim.fn, nsim= 1000, .progress="win")


## back-transformation needed!
# estim.GesamtscoreLog.df <- cbind (estim.df, exp (extract.ci (GesamtscoreLog.boot)))
# estim.val <- exp(estim.val)-0.5   #Zähldaten /Log Transformation
estim.GesamtscoreLog.df <- cbind (estim.df, exp(extract.ci (GesamtscoreLog.boot))-0.5)
summary(estim.GesamtscoreLog.df)

estim.GesamtscoreLog.df

### Graphical presentation ------------
boxplot(Score ~ Status + Phase, health.df)

# ADWG Model --------------------------------------------------
ADWG.full <- blmer (ADWG ~
                       ## main effects
                       Status1 +
                       Phase1 + Phase2 + Phase3 + Phase4 +
                       Phase5 + 
                       ## Two-way interaction
                       Status1:Phase1 + Status1:Phase2 +
                       Status1:Phase3 + Status1:Phase4 +
                       Status1:Phase5 + 
                       ## Random effect
                       (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                       data = weight.df,  REML = FALSE)
summary (ADWG.full)

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

#Show outliers
weight.df[resid (ADWG.full)>  1.0,]
weight.df[resid (ADWG.full)< -(1.5),]
weight.df[fitted.values (ADWG.full)< -2.0,]

plotResiduals (ADWG.resid, weight.df[, 'Status'])
boxplot (split (resid (ADWG.full), weight.df [, 'Status']))

plotResiduals (ADWG.resid, weight.df[, 'Phase'])            
boxplot (split (resid (ADWG.full), weight.df [, 'Phase']))

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

## zero model
ADWG.zero <- blmer (ADWG ~ 1 +
                    (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                    data = weight.df, REML = FALSE)
summary(ADWG.zero)

## single effects
ADWG.status <- blmer (ADWG ~
                         ## main effects
                         #                        Status1 +
                         Phase1 + Phase2 + Phase3 + Phase4 +
                         Phase5 + 
                         ## Two-way interaction
                         Status1:Phase1 + Status1:Phase2 +
                         Status1:Phase3 + Status1:Phase4 +
                         Status1:Phase5 + 
                         ## Random effect
                         (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                         data = weight.df, REML = FALSE)

ADWG.phase <- blmer (ADWG ~
                        ## main effects
                        Status1 +
                        #                           Phase1 + Phase2 + Phase3 + Phase4 + ...
                        ## Two-way interaction
                        Status1:Phase1 + Status1:Phase2 +
                        Status1:Phase3 + Status1:Phase4 +
                        Status1:Phase5 + 
                        ## Random effect
                        (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                        data = weight.df, REML = FALSE)

ADWG.int <- blmer (ADWG ~
                      ## main effects
                      Status1 +
                      Phase1 + Phase2 + Phase3 + Phase4 +
                      Phase5 + 
                      ## Two-way interaction
                      #                         Status1:Phase1 + Status1:Phase2 +
                      #                         Status1:Phase3 + Status1:Phase4 +
                      ## Random effect
                      (1 | LargeGroup/Batch/FosterCow/CalfID/Phase) + (1 | location) + (1 | CalWeek),
                      data = weight.df, REML = FALSE)

### ANOVA model comparison (quick but anti-conservative) ---------
# global
anova (ADWG.full, ADWG.zero)

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

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

save.image()


### Bootstrap ADWG --------------------------------------------------------
# global
PBmodcomp.ADWG.zero <- PBmodcomp (ADWG.full, ADWG.zero)

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

PBmodcomp.ADWG.phase <- PBmodcomp(ADWG.full, ADWG.phase)

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

#Summaries
summary (PBmodcomp.ADWG.zero)
summary (PBmodcomp.ADWG.status)
summary (PBmodcomp.ADWG.phase)
summary (PBmodcomp.ADWG.int)


### Model estimate --------------
estim.df <- data.frame (Status1= rep (c (1, -1), 6),
                        Phase1= rep (c (1, 0, 0, 0, 0, -1), rep(2,6)),
                        Phase2= rep (c (0, 1, 0, 0, 0, -1), rep(2,6)),
                        Phase3= rep (c (0, 0, 1, 0, 0, -1), rep(2,6)),
                        Phase4= rep (c (0, 0, 0, 1, 0, -1), rep(2,6)),
                        Phase5= rep (c (0, 0, 0, 0, 1, -1), rep(2,6)),
                        Phase6= rep (c (0, 0, 0, 0, 0, -1), rep(2,6)),
                        Sex1 =0)
## estimate and confidence interval bootstrap
estim.fn <- function (x) predict (x, newdat= estim.df, re.form= NA)
ADWG.boot <- bootMer (ADWG.full, estim.fn, nsim= 1000, .progress="win")

 
## back-transformation ? not needed here
# estim.ADWG.df <- cbind (estim.df, exp (extract.ci (ADWG.boot)))
estim.ADWG.df <- cbind (estim.df, extract.ci (ADWG.boot))
summary(estim.ADWG.df)

estim.ADWG.df


#### Status Estimate ----------
## Model estimate
estim_status.df <- data.frame (Status1= c (1, -1),
                               Phase1= 0,
                               Phase2= 0,
                               Phase3= 0,
                               Phase4= 0,
                               Phase5= 0,
                               Phase6= 0,
                               Sex1 =0)
## estimate and confidence interval bootstrap
estim.fn <- function (x) predict (x, newdat= estim_status.df, re.form= NA)
ADWG.boot_status <- bootMer (ADWG.full, estim.fn, nsim= 1000, .progress="win")


## back-transformation ? not needed here
# estim.ADWG.df <- cbind (estim.df, exp (extract.ci (ADWG.boot)))
estim.ADWG_status.df <- cbind (estim.df, extract.ci (ADWG.boot_status))
summary(estim.ADWG_status.df)
estim.ADWG_status.df
### Graphical presentation----------------
boxplot (ADWG ~ Status + Phase, weight.df)

# Final Weight Models --------------------------------------------------
##Model Final weight
final.weight.df$finalweight <-final.weight.df$Weight

finalweight.full <- blmer (log(finalweight) ~
                      ## main effects
                      Status1 +
                      Age_tf +
                      (1 | LargeGroup/Batch/FosterCow),
                    data = final.weight.df,
                    #control = lmerControl(optimizer="bobyqa"),
                    REML = FALSE)
summary (finalweight.full)

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

#Show outliers
#final.weight.df[resid (finalweight.full)>  1.5,]
#final.weight.df[resid (finalweight.full)< -1.6,]
#final.weight.df[fitted.values (finalweight.full)< -2.0,]

plotResiduals (finalweight.resid, final.weight.df[, 'Status'])
boxplot (split (resid (finalweight.full), final.weight.df [, 'Status']))

plotResiduals (finalweight.resid, final.weight.df[, 'Age_tf'])
boxplot (split (resid (finalweight.full), final.weight.df [, 'Age_tf']))

qqnorm (ranef (finalweight.full) [['Batch:LargeGroup']] [, 1])
qqnorm (ranef (finalweight.full) [['FosterCow:(Batch:LargeGroup)']] [, 1])

## zero model
finalweight.zero <- blmer (log(finalweight) ~ 1 +
                      (1 | LargeGroup/Batch/FosterCow),
                    data = final.weight.df, 
                    #control = lmerControl(optimizer="bobyqa"),
                    REML = FALSE)
summary(finalweight.zero)


## single effects
finalweight.status <- blmer (log(finalweight) ~
                        ## main effects
                        #                        Status1 +
                        Age_tf +
                        ## Random effect
                        (1 | LargeGroup/Batch/FosterCow),
                        data = final.weight.df, 
                        #control = lmerControl(optimizer="bobyqa"),
                        REML = FALSE)
summary(finalweight.status)

finalweight.age <- blmer (log(finalweight) ~
                               ## main effects
                               Status1 +
                               #Age_tf +
                               ## Random effect
                               (1 | LargeGroup/Batch/FosterCow),
                             data = final.weight.df, 
                             #control = lmerControl(optimizer="bobyqa"),
                          REML = FALSE)
summary(finalweight.age)

### Anovas Final Weight -----------------------------------------------------
## model comparison (quick but anti-conservative)
# global
anova (finalweight.full, finalweight.zero)

# main effects
anova (finalweight.full, finalweight.status)
anova (finalweight.full, finalweight.age)

save.image()


### Bootstrap finalweight --------------------------------------------------------
# global
PBmodcomp.finalweight.zero <- PBmodcomp (finalweight.full, finalweight.zero)

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

PBmodcomp.finalweight.age <- PBmodcomp (finalweight.full, finalweight.age)

#Summaries
summary (PBmodcomp.finalweight.zero)
summary (PBmodcomp.finalweight.status)
summary (PBmodcomp.finalweight.age)

### Model estimate------------------
estim.df <- data.frame (Status1= c (1, -1),
                        Age_tf= 0,
                        Sex1=0)

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

## back-transformation ? yes!
# estim.finalweight.df <- cbind (estim.df, exp (extract.ci (finalweight.boot)))
estim.finalweight.df <- cbind (estim.df, exp(extract.ci (finalweight.boot)))
summary(estim.finalweight.df)
estim.finalweight.df

### Graphical presentation------------------------
boxplot (Weight ~ Status, final.weight.df)

# Figures GERMAN ---------------------------------------------------------
###Gesundheit Boxplot---------
#cols <- rep(c("#4F81BD", "#9BCB59"), times = 7)  # AK blue, EK green
cols <- rep(c("#1F4E79", "#A9D18E"), times = 7)  # AK blue, EK green - more contrast -> print

jpeg("Gesundheit Boxplot.jpg", width = 18, height = 13, units = "cm", res=300, pointsize=11)
par(bty = 'n', las = 1
    #, mfrow = c (4, 4)
    , mar = c(7, 4.5, 0.5, 0.5)
    #, oma = c(6, 6, 0, 0))
)
boxplot (Score ~ Status + Phase, health.df,
         ylab= '', xlab= '', yaxt= 'n', xaxt= 'n',
         whisklty= 1, range= 0, ylim=c(0,8), col=cols)
axis (2, seq(0,8,1), seq(0,8,1))
title(ylab = "Erkrankungsscore", cex.lab = 1.2)
for (i in seq (1, 13, 2)) axis (1, c (i, i+1), c ('AK', 'EK'))
axis (1, seq (1.5, 13.5, 2), c ('Klein-\ngruppe', '','Großgruppe\n','',
                                'Absetzen\n(1)', 'Absetzen\n(2)', 'Abgesetzt\n'),
      tick= FALSE, line= 3.5)
axis (1, seq(1.5, 13.5, 2), c ('','(Anfang)', '(Mitte)',
                               '(Ende)','','',''),
      tick= FALSE, line= 3.5)
for (i in seq (1, 13, 2)) {
  lines (c (i, i+1), estim.GesamtscoreLog.df [c (i, i+1), 'estim'], lwd= 4)
  lines (c (i, i+1), estim.GesamtscoreLog.df [c (i, i+1), 'lo.ci'], lwd= 2)
  lines (c (i, i+1), estim.GesamtscoreLog.df [c (i, i+1), 'up.ci'], lwd= 2)
}
rm (i)
dev.off()

estim.GesamtscoreLog.df

### Tageszunahmen------------------
jpeg("Tageszunahmen_ohneStrich.jpg", width = 18, height = 13, units = "cm", res=300, pointsize = 11)
par(bty = 'n', las = 1
    #, mfrow = c (4, 4)
    , mar = c(6, 5, 0.5, 0.5)
    #, oma = c(6, 6, 0, 0))
)
boxplot (ADWG ~ Status + Phase, weight.df,
         ylab= '', xlab= '', yaxt= 'n', xaxt= 'n',
         whisklty= 1, range= 0,ylim = c(-1, 2.5), col=cols)
#abline(h = seq(-2,2.5,0.5), col = "lightgrey", lty="dashed")
#abline(h =0.85,col= "darkred", lty=1, cex=10)
title(ylab = "Tageszunahmen [g]", cex.lab = 1.2, line=3.5)
axis (2, seq(-1,2.5,0.5), seq(-1000,2500,500))
for (i in seq (1, 13, 2)) axis (1, c (i, i+1), c ('AK', 'EK'))
axis (1, seq (1.5, 13.5, 2), c ('Klein-\ngruppe', '','Großgruppe\n','',
                                'Absetzen\n(1)', 'Absetzen\n(2)', 'Abgesetzt\n'),
      tick= FALSE, line= 3.5)
axis (1, seq(1.5, 13.5, 2), c ('','(Anfang)', '(Mitte)',
                               '(Ende)','','',''),
      tick= FALSE, line= 3.5)

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

estim.ADWG.df

### Absetzgewicht------------
jpeg("AbsetzWeight.jpg", width = 16, height = 13, units = "cm", res=300, pointsize = 11)
par(bty = 'n', las = 1
    #, mfrow = c (4, 4)
    , mar = c(4, 5, 1, 0.5)
    #, oma = c(6, 6, 0, 0))
)
boxplot (Weight ~ Status, final.weight.df,
         ylab= '', xlab= '', yaxt= 'n', xaxt= 'n',
         whisklty= 1, range= 0,ylim = c(100, 250), col=cols)
title(ylab = "AbsetzWeight [kg]", cex.lab = 1.2, line=3.5)
abline(h = seq(100, 250,25), col = "grey30", lty="dashed")
axis (2, seq(100, 250,50))
for (i in seq (1, 11, 2)) axis (1, c (i, i+1), c ('Ammenkälber', 'Eigene Kälber'))
for (i in seq (1, 11, 2)) {
  lines (c (i, i+1), estim.finalweight.df [c (i, i+1), 'estim'], lwd= 4)
  lines (c (i, i+1), estim.finalweight.df [c (i, i+1), 'lo.ci'], lwd= 2)
  lines (c (i, i+1), estim.finalweight.df [c (i, i+1), 'up.ci'], lwd= 2)
}
rm (i)
dev.off()
estim.finalweight.df

#Figures ENGLISH JDS--------
## Graphical presentation
#cols <- rep(c("#4F81BD", "#9BCB59"), times = 7)  # AK blue, EK green
cols <- rep(c("#1F4E79", "#A9D18E"), times = 7)  # AK blue, EK green - more contrast -> print
col_bw <- rep (c("gray30", "gray80"))


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

# Sample size of Calf-Pairs per Phase
n.pairs_w <- tapply(weight.df$FosterCow,
                    weight.df$Phase,
                    function(x) length(unique(na.omit(x))))

### Health Boxplot---------
jpeg("Health Score Boxplot_english_JDS.jpg", width = 14, height = 10, units = "cm", res=300, pointsize=8)
par(bty = 'n', las = 1
    #, mfrow = c (4, 4)
    , mar = c(7, 4.5, 0.5, 0.5)
    #, oma = c(6, 6, 0, 0))
)
boxplot (Score ~ Status + Phase, health.df,
         ylab= '', xlab= '', yaxt= 'n', xaxt= 'n',
         whisklty= 1, range= 0, ylim=c(-0.5,9), col=cols)
axis (2, seq(0,9,1), seq(0,9,1))
title(ylab = "Health Impairment Score", cex.lab=1.2)
for (i in seq (1, 13, 2)) axis (1, c (i, i+1), c ('FC', 'OC'))
axis (1, seq (1.5, 13.5, 2), c ('small group\n\n', '','large group\n\n','',
                                'weaning\n(1)\n', 'weaning\n(2)\n', 'weaned\n\n'),
      tick= FALSE, line= 3.5)
axis (1, seq(1.5, 13.5, 2), c ('','(beginning)', '(mid)',
                               '(end)','','',''),
      tick= FALSE, line= 3.0)

axis(1, seq(1.5, 13.5, 2),
     paste0('(n = ', n.pairs, ' pairs)'),
     tick = FALSE, line = -2.25, cex.axis = 0.75)

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

### Weight gain Boxplot------------------
jpeg("ADWG_english_JDS.jpg", width = 14, height = 10, units = "cm", res=300, pointsize = 8)
par(bty = 'n', las = 1
    #, mfrow = c (4, 4)
    , mar = c(6, 5, 0.5, 0.5)
    #, oma = c(6, 6, 0, 0))
)
boxplot (ADWG ~ Status + Phase, weight.df,
         ylab= '', xlab= '', yaxt= 'n', xaxt= 'n',
         whisklty= 1, range= 0,ylim = c(-0.7, 2.5), col=cols)
#abline(h = seq(-2,2.5,0.5), col = "lightgrey", lty="dashed")
#abline(h =0.85,col= "darkred", lty=1, cex=10)
title(ylab = "ADWG [g]", cex.lab = 1.2, line=3.5)
axis (2, seq(-0.5,2.5,0.5), seq(-500,2500,500))
for (i in seq (1, 13, 2)) axis (1, c (i, i+1), c ('FC', 'OC'))
axis (1, seq (1.5, 13.5, 2), c ('small group\n\n', '','large group\n\n','',
                                'weaning\n(1)\n', 'weaning\n(2)\n', 'weaned\n\n'),
      tick= FALSE, line= 3.5)
axis (1, seq(1.5, 13.5, 2), c ('','(beginning)', '(mid)',
                               '(end)','','',''),
      tick= FALSE, line= 3.0)
axis(1, seq(1.5, 11.5, 2),
     paste0('(n = ', n.pairs_w, ' pairs)'),
     tick = FALSE, line = -2.25, cex.axis = 0.75)

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

### Weaning (final) Weight boxplot ------------
jpeg("final weight_english_JDS.jpg", width = 8.9, height = 10, units = "cm", res=300, pointsize = 8)
par(bty = 'n', las = 1
    #, mfrow = c (4, 4)
    , mar = c(4, 5, 1, 0.5)
    #, oma = c(6, 6, 0, 0))
)
boxplot (Weight ~ Status, final.weight.df,
         ylab= '', xlab= '', yaxt= 'n', xaxt= 'n',
         whisklty= 1, range= 0,ylim = c(100, 250), col=cols)
title(ylab = "Weaning Weight [kg]", cex.lab = 1.2, line=3.5)
#abline(h = seq(100, 250,25), col = "grey30", lty="dashed")
axis (2, seq(100, 250,50))
for (i in seq (1, 2, 2)) axis (1, c (i, i+1), c ('Foster Calves', 'Own Calves'))
for (i in seq (1, 2, 2)) {
  lines (c (i, i+1), estim.finalweight.df [c (i, i+1), 'estim'], lwd= 4)
  lines (c (i, i+1), estim.finalweight.df [c (i, i+1), 'lo.ci'], lwd= 2)
  lines (c (i, i+1), estim.finalweight.df [c (i, i+1), 'up.ci'], lwd= 2)
}
rm (i)
dev.off()
estim.finalweight.df

##XXX Additional XXX---------------------
###   Prevalences -----------------------
nrow(health.df[health.df$Vitality==1,])
nrow(health.df[health.df$Vitality==1,])/nrow(health.df)
summary(health.df[health.df$Vitality==1,])

health.df$Resp <- as.numeric(health.df$Coughing==1 |
                               health.df$Breathing==1 |
                               health.df$Nasal_Discharge==1 |
                               health.df$Ocular_Discharge==1)

health.df$Vit_Resp <- as.numeric(health.df$Vitality==1 &
                                    (health.df$Coughing==1 |
                                       health.df$Breathing==1 |
                                       health.df$Nasal_Discharge==1 |
                                       health.df$Ocular_Discharge==1))
health.df$Vit_Diarr <- as.numeric(health.df$Vitality==1 &
                                     health.df$Diarrhea==1 )

summary(health.df$Vit_Diarr)
summary(health.df$Vit_Resp)

#Health impairment Score >2 
nrow(health.df)
nrow(health.df[health.df$Score>=2,])
nrow(health.df[health.df$Score>=2,])/nrow(health.df)

#Diarrhea + Allg
nrow(health.df[health.df$Diarrhea == 1, ])
nrow(health.df[health.df$Diarrhea == 1, ])/nrow(health.df)

nrow(health.df[health.df$Diarrhea == 1 &
               health.df$Vitality == 1, ])

nrow(health.df[health.df$Diarrhea == 1 &
                 health.df$Vitality == 1, ])/nrow(health.df[health.df$Diarrhea == 1, ])

#Respiratory + Vitality
nrow(health.df[health.df$Resp == 1, ])
nrow(health.df[health.df$Resp == 1, ])/nrow(health.df)

nrow(health.df[health.df$Resp == 1 &
                 health.df$Vitality == 1, ])

nrow(health.df[health.df$Resp == 1 &
                 health.df$Vitality == 1, ])/nrow(health.df[health.df$Resp == 1, ])


summary_prevalences <- aggregate(cbind(Vitality,
                                               BCS,Coat_Condition,Diarrhea,Coughing,Nasal_Discharge,
                                               Ocular_Discharge,Breathing, Resp, Vit_Resp,Vit_Diarr) 
                                         ~ Status + Phase, 
                                         data = health.df, 
                                         FUN = function(x) round(mean(x == 1, na.rm = TRUE)*100, digits = 2))
summary(summary_prevalences)

write.csv(summary_prevalences, file = "Prevalences_Healthscore_Status_and_Phase.csv", row.names = FALSE)

#### Plots -------------------
cols <- rep(c("#1F4E79", "#A9D18E"), times = 7)  # AK blue, EK green - more contrast -> print

#Impaired Vitality + Respiratory symptoms
par(bty = 'n', las = 1,
    mar = c(5, 6, 1, 1)) #, mfrow = c (4, 4), , oma = c(6, 8, 0, 0))
barplot (Vit_Resp ~ Status + Phase, summary_prevalences,
         ylab= 'Prevalence of impaired vitality\nwith respiratory symptoms (%)',
         ylim=c(0,30), col=cols, beside=TRUE, width=1)
axis (2, seq(0,30,10), seq(0,30,10))
legend("topright",
       legend = levels(summary_prevalences$Status),
       pch= 15,
       col = cols,
       bty = "o",   # Kasten
       cex = 0.8,
       pt.cex = 1.2,
       y.intersp = 1.2)

#Impaired Vitality + Diarrhea symptoms
par(bty = 'n', las = 1,
    mar = c(5, 6, 1, 1)) #, mfrow = c (4, 4), , oma = c(6, 8, 0, 0))
barplot (Vit_Diarr ~ Status + Phase, summary_prevalences,
         ylab= 'Prevalence of impaired\nvitality with diahrrea (%)',
         ylim=c(0,30), col=cols, beside=TRUE, width=1)
axis (2, seq(0,30,10), seq(0,30,10))
legend("topright",
       legend = levels(summary_prevalences$Status),
       pch= 15,
       col = cols,
       bty = "o",   # Kasten
       cex = 0.8,
       pt.cex = 1.2,
       y.intersp = 1.2)


summary_prevalences_long <- summary_prevalences %>%
  pivot_longer(
    cols = -c(Status, Phase),  # All rows except Status and Phase converted
    names_to = "Parameter",    # Name of new column for variables
    values_to = "Prevalence"   # Name of new column for prevalence-values
  )

#Prevalences per Indicator
ggplot(summary_prevalences_long, aes(x = Phase, y = Prevalence, fill=Status)) +
  geom_col(position = "dodge") +  
  facet_wrap(~ Parameter, scales = "free_y") +  
  labs(
    title = "Prevalence of health indicators according to Status and Phase",
    x = "Health indicator",
    y = "Prevalence (%)",
    fill = "Status"
  ) +
  scale_fill_manual(
    values = c("foster calf" = "#1F4E79", "own calf" = "#A9D18E"),
    name = "Status"
  ) +
  scale_y_continuous(
    limits = c(0, 100),                    
    breaks = seq(0, 100, by = 25)) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),  
    legend.position = "bottom"
  )

# XXX Quit Session in R XXX -------------------------------------------------------------
save.image ()
q ('no')