09_Figure2_3_Best_models.R

####################################################################################################################
# Analysis of implant data from foxes for:
# Shapiro et al.
# Potential for real-time health and welfare monitoring in experimental rabies infection in red fox (Vulpes vulpes) 
# using implants
# 
# Script 9: Script to create Figures 2 and 3 illustrating best model for each fox 
#          
#
# Uses data frames created for each fox in scripts 05 - 08 :
#    avant_0fae.act.24h, apres_0fae.act.24h, avant_0faf.act.24h, apres_0faf.act.24h,
#    avant_0fb1.act.24h, apres_0fb1.act.24h, avant_0fb0.act.24h, apres_0fb0.act.24h 
# 
# Script tested for R version 4.1.2
####################################################################################################################


####################################################################################################################
# Load packages ####
####################################################################################################################
# Load packages (install if necessary)
library(tidyverse)
library(qcc)
library(ggpubr)
####################################################################################################################

####################################################################################################################
# Function : Reruns best model with best parameter settings ####
####################################################################################################################

# Define number of sigmas, must be done for each fox (see below) :
nsigma=

# Function :
data_prep_ewma_ggplot <- function(data_calibration,
                                  data_test){
  # Create variables
  n_calib <- nrow(data_calibration)
  n_total <- nrow(data_calibration) + nrow(data_test)
  
  # Create QCC object
  # MAKE SURE TO CHANGE THE NSIGMAS FOR EACH FOX
  q <- ewma(data_calibration, lambda = 0.3, nsigmas = nsigma, 
                                  newdata = data_test, 
                                 plot=F)
  
  # Create dataframe with needed information
  qcc_data <- data.frame(
    Date = seq(1,(length(q$statistics) + length(q$newstats)), 1),
    Value = rep(NA, n_total), #vector of NA, same length as total observations
    UCL = q$limits[, 2],      #get UCL from QCC
    LCL = q$limits[, 1],
    EWMA = q$y)      #get LCL from QCC
    
  #Nb_values = rep(NA, n_total))  #vector of NA, same length as total observations
  
  # Adapt dataframe
  qcc_data$Value[1:n_calib] <- q$statistics     #fill Value variable
  qcc_data$Value[(n_calib+1):n_total] = q$newstats
  qcc_data$Date[1:n_calib] <- rownames(q$data)  #fill Date variable
  qcc_data$Date[(n_calib+1):n_total] <- rownames(q$newdata)
  qcc_data$Date <- factor(qcc_data$Date, levels = qcc_data$Date) #Date as factor
  qcc_data$EWMA = q$y
  
  # Make a TRUE/FALSE variable if value is outside the limits
  qcc_data$alert <- with(qcc_data, EWMA > UCL | EWMA < LCL)
 
  return(qcc_data)
  
}
####################################################################################################################


####################################################################################################################
# Data preparation for figure ####
####################################################################################################################

# Create sequence of dates for missing values to add to figures
missingdates.bestmod <- seq(as.Date('2022-07-28'),as.Date('2022-08-08'),by = 'day') %>%
  as.data.frame() %>%
  rename(., Date=.)
####################################################################################################################


####################################################################################################################
# Apply function to each fox ####
####################################################################################################################


####################################################################################################################
# Fox 0fae ####
####################################################################################################################

# Apply the data prep function above
# First define nsigmas for model
nsigma = 2  

# Apply function and create table from results to make figure
best.0fae.mod <- data_prep_ewma_ggplot(avant_0fae.act.24h, apres_0fae.act.24h) %>%
  mutate(alert=ifelse(alert==TRUE, 'Alarm','Normal range'),
         Date=as.Date(Date)) %>%
  # Add missing days
  bind_rows(missingdates.bestmod)

# Assign rownames
rownames(best.0fae.mod) <- best.0fae.mod$Date

# Create the figure
best.0fae.mod.plot <- ggplot(best.0fae.mod, aes(x=Date)) +
  #ylim(0,2.5) +
  # Lines
  # Horizontal line for mean
  geom_hline(aes(yintercept = ewma2s.0fae.act.24h$center), 
             color = "black", linetype = "solid") + 
  # Line for EWMA values
  geom_line(aes(y = EWMA, group = 1), color = "black") +
  # Upper limit
  geom_line(aes(y = UCL, group = 1), color = "black", linetype="dashed", size=1.2) + #UCL
  # Lower limite
  geom_line(aes(y = LCL, group = 1), color = "black", linetype="dashed", size=1.2) + #LCL
  # Vertical line for end of asymptomatic / calibration period
  geom_vline(xintercept = as.Date('2022-07-19'), linetype="dotdash", 
             color = "blue4", size=1) +
  # Vertical lines for days with missing data
  geom_vline(xintercept = as.Date('2022-07-28'), linetype="dotted", 
             color = "azure4", size=1.8) +
  geom_vline(xintercept = as.Date('2022-08-08'), linetype="dotted", 
             color = "azure4", size=1.8) +
  # Vertical line for day of death
  geom_vline(xintercept = as.Date('2022-08-18'), linetype="longdash",
             color = "red4", size=0.5) +
  # Dots
  # Dots for EWMA statistic points, color by alarm or not
  geom_point(aes(y = EWMA, color = alert), size = 2.5) +
  # Crosses for actual mean value
  geom_point(aes(y = Value), shape=3, color="black",size=4) +
  
  # Legend 
  scale_color_manual(values = c(
    "Alarm" = "red",
    "Normal range" = "black"),na.translate = F) +
  
  # Labels
  labs(
    x = "Date",
    y = 'Activity (g)') +
  
  #color = c("") +
  theme_classic() +
  scale_x_datetime(
    date_breaks = '2 days',
    # use anonymous function instead of character vector (works for all scales)
    label = \(datetime) ymd(datetime)) +
  theme(legend.text=element_text(size=17),
        text = element_text(size = 17), axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(legend.title=element_blank()) 
####################################################################################################################


####################################################################################################################
# For 0faf ####
####################################################################################################################

# Apply the data prep function above
# First define nsigmas for model
nsigma = 3  

# Apply function and create table from results to make figure
best.0faf.mod <- data_prep_ewma_ggplot(avant_0faf.act.24h, apres_0faf.act.24h) %>%
  mutate(alert=ifelse(alert==TRUE, 'Alarm','Normal range'),
         Date=as.Date(Date)) %>%
  # Add missing days
  bind_rows(missingdates.bestmod)

# Rownames (probably not necessary)
rownames(best.0faf.mod) <- best.0faf.mod$Date


# Make the figure
best.0faf.mod.plot <- ggplot(best.0faf.mod, aes(x=Date)) +
  #ylim(0,2.5) +
  # Lines
  # Horizontal line for mean
  geom_hline(aes(yintercept = ewma3s.0faf.act.24h$center), 
             color = "black", linetype = "solid") + 
  # Line for EWMA values
  geom_line(aes(y = EWMA, group = 1), color = "black") +
  # Upper limit
  geom_line(aes(y = UCL, group = 1), color = "black", linetype="dashed", size=1.2) + #UCL
  # Lower limite
  geom_line(aes(y = LCL, group = 1), color = "black", linetype="dashed", size=1.2) + #LCL
  # Vertical line for end of asymptomatic / calibration period
  geom_vline(xintercept = as.Date('2022-07-19'), linetype="dotdash", 
             color = "blue4", size=1) +
  # Vertical line for day of death
  geom_vline(xintercept = as.Date('2022-08-08'), linetype="longdash",
             color = "red4", size=0.5) +
  # Vertical line for missing data
  geom_vline(xintercept = as.Date('2022-07-28'), linetype="dotted", 
             color = "azure4", size=1.8) +
  geom_vline(xintercept = as.Date('2022-08-08'), linetype="dotted", 
             color = "azure4", size=1.8) +
  # Dots
  # Dots for EWMA statistic points, color by alarm or not
  geom_point(aes(y = EWMA, color = alert), size = 2.5) +
  # Crosses for actual mean value
  geom_point(aes(y = Value), shape=3, color="black",size=4) +
  
  # Legend 
  scale_color_manual(values = c(
    "Alarm" = "red",
    "Normal range" = "black"),na.translate = F) +
  
  # Labels
  labs(
    x = "Date",
    y = 'Activity (g)') +
  
  #color = c("") +
  theme_classic() +
  scale_x_datetime(
    date_breaks = '2 days',
    # use anonymous function instead of character vector (works for all sclaes)
    label = \(datetime) ymd(datetime)) +
  #scale_x_datetime(date_breaks = "2 days", date_labels = "%Y/%m/%d") +
  theme(legend.text=element_text(size=17),
        text = element_text(size = 17), axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(legend.title=element_blank())
####################################################################################################################


####################################################################################################################
# For 0fb1 ####
####################################################################################################################

# Apply the data prep function above
# First define nsigmas for model
nsigma = 4

# Apply function and create table from results to make figure
best.0fb1.mod <- data_prep_ewma_ggplot(night.avant.0fb1.act, night.apres.0fb1.act) %>%
  mutate(alert=ifelse(alert==TRUE, 'Alarm','Normal range'),
         Date=as.Date(Date)) %>%
  # Add missing days
  bind_rows(missingdates.bestmod) %>%
  filter(Date <= "2022-08-05")

# Assign rownames
rownames(best.0fb1.mod) <- best.0fb1.mod$Date


# Make the figure
best.0fb1.mod.plot <- ggplot(best.0fb1.mod, aes(x=Date)) +
  #ylim(0,2.5) +
  # Lines
  # Horizontal line for mean
  geom_hline(aes(yintercept = ewma4s.0fb1.act.night$center), 
             color = "black", linetype = "solid") + 
  # Line for EWMA values
  geom_line(aes(y = EWMA, group = 1), color = "black") +
  # Upper limit
  geom_line(aes(y = UCL, group = 1), color = "black", linetype="dashed", size=1.2) + #UCL
  # Lower limite
  geom_line(aes(y = LCL, group = 1), color = "black", linetype="dashed", size=1.2) + #LCL
  # Vertical line for end of asymptomatic / calibration period
  geom_vline(xintercept = as.Date('2022-07-19'), linetype="dotdash", 
             color = "blue4", size=1) +
  # Vertical line for day of death
  geom_vline(xintercept = as.Date('2022-08-04'), linetype="longdash",
             color = "red4", size=0.5) +
  # Dotted vertical lines for data gap
  geom_vline(xintercept = as.Date('2022-07-27'), linetype="dotted", 
             color = "azure4", size=1.8) +
  geom_vline(xintercept = as.Date('2022-08-04'), linetype="dotted", 
             color = "azure4", size=1.8) +
  # Dots
  # Dots for EWMA statistic points, color by alarm or not
  geom_point(aes(y = EWMA, color = alert), size = 2.5) +
  # Crosses for actual mean value
  geom_point(aes(y = Value), shape=3, color="black",size=4) +
  
  # Legend 
  scale_color_manual(values = c(
    "Alarm" = "red",
    "Normal range" = "black"),na.translate = F) +
  
  # Labels
  labs(
    x = "Date",
    y = 'Activity (g)') +
  
  #color = c("") +
  theme_classic() +
  #scale_x_datetime(date_breaks = "2 days", date_labels = "%Y/%m/%d") +
  scale_x_datetime(
    date_breaks = '2 days',
    # use anonymous function instead of character vector (works for all sclaes)
    label = \(datetime) ymd(datetime)) +
  theme(legend.text=element_text(size=17),
        text = element_text(size = 17), axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(legend.title=element_blank())
####################################################################################################################


####################################################################################################################
# 0fb0 ####
####################################################################################################################

# Apply the data prep function above
# First define nsigmas for model
nsigma = 4

# Apply function and create table from results to make figure
best.0fb0.mod <- data_prep_ewma_ggplot(morning.avant.0fb0.act, morning.apres.0fb0.act) %>%
  mutate(alert=ifelse(alert==TRUE, 'Alarm','Normal range'),
         Date=as.Date(Date)) %>%
  # Add missing days
  bind_rows(missingdates.bestmod) %>%
  filter(Date <='2022-08-19') %>%
  dplyr::arrange(Date)

# Assign rownames
rownames(best.0fb0.mod) <- best.0fb0.mod$Date 


# Make first figure for Figure 2 up to August 18 
# Lines
# Horizontal line for mean
best.0fb0.mod.plot <- ggplot(best.0fb0.mod, aes(x=Date)) +
  #ylim(0,2.5) +
  # Lines
  # Horizontal line for mean
  geom_hline(aes(yintercept = ewma4s.0fb0.act.morning$center), 
             color = "black", linetype = "solid") + 
  # Line for EWMA values
  geom_line(aes(y = EWMA, group = 1), color = "black") +
  # Upper limit
  geom_line(aes(y = UCL, group = 1), color = "black", linetype="dashed", size=1.2) + #UCL
  # Lower limite
  geom_line(aes(y = LCL, group = 1), color = "black", linetype="dashed", size=1.2) + #LCL
  # Vertical line for end of asymptomatic / calibration period
  geom_vline(xintercept = as.Date('2022-07-19'), linetype="dotdash", 
             color = "blue4", size=1) +
  geom_vline(xintercept = as.Date('2022-07-28'), linetype="dotted", 
             color = "azure4", size=1.8) +
  geom_vline(xintercept = as.Date('2022-08-08'), linetype="dotted", 
             color = "azure4", size=1.8) +
  # Dots
  # Dots for EWMA statistic points, color by alarm or not
  geom_point(aes(y = EWMA, color = alert), size = 2) +
  # Crosses for actual mean value
  geom_point(aes(y = Value), shape=3, color="black",size=3) +
  
  # Legend 
  scale_color_manual(values = c(
    "Alarm" = "red",
    "Normal range" = "black"),na.translate = F) +
  
  # Labels
  labs(
    x = "Date",
    y = 'Activity (g)') +
  
  theme_classic() +
  # This was the only way to get the correct labels for the dates
  scale_x_datetime(
    date_breaks = '2 days',
    # use anonymous function instead of character vector (works for all sclaes)
    label = \(datetime) ymd(datetime)) +
  theme(legend.text=element_text(size=17),
        text = element_text(size = 17), axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(legend.title=element_blank()) +
  # This needed to get the right range of dates for the labels
  expand_limits(x =c(min(best.0fb0.mod$Date)+1, max(best.0fb0.mod$Date)))


# Figure 3 for 0fb0, all data
best.0fb0.mod.full <- data_prep_ewma_ggplot(morning.avant.0fb0.act, morning.apres.0fb0.act) %>%
  mutate(alert=ifelse(alert==TRUE, 'Alarm','Normal range'),
         Date=as.Date(Date)) %>%
  # Add missing days
  bind_rows(missingdates.bestmod)

best.0fb0.mod.full.plot <- ggplot(best.0fb0.mod.full, aes(x=Date)) +
  # Lines
  # Horizontal line for mean
  geom_hline(aes(yintercept = ewma4s.0fb0.act.morning$center), 
             color = "black", linetype = "solid") + 
  # Line for EWMA values
  geom_line(aes(y = EWMA, group = 1), color = "black") +
  # Upper limit
  geom_line(aes(y = UCL, group = 1), color = "black", linetype="dashed", size=1.2) + #UCL
  # Lower limite
  geom_line(aes(y = LCL, group = 1), color = "black", linetype="dashed", size=1.2) + #LCL
  # Vertical line for end of asymptomatic / calibration period
  geom_vline(xintercept = as.Date('2022-07-19'), linetype="dotdash", 
             color = "blue4", size=1) +
  geom_vline(xintercept = as.Date('2022-07-28'), linetype="dotted", 
             color = "azure4", size=1.8) +
  geom_vline(xintercept = as.Date('2022-08-08'), linetype="dotted", 
             color = "azure4", size=1.8) +
  # Dots
  # Dots for EWMA statistic points, color by alarm or not
  geom_point(aes(y = EWMA, color = alert), size = 2) +
  # Crosses for actual mean value
  geom_point(aes(y = Value), shape=3, color="black",size=3) +
  
  # Legend 
  scale_color_manual(values = c(
    "Alarm" = "red",
    "Normal range" = "black"),na.translate = F) +
  # Labels
  labs(
    x = "Date",
    y = 'Activity (g)') +
  
  #color = c("") +
  theme_classic() +
  scale_x_datetime(
    date_breaks = '2 days',
    # use anonymous function instead of character vector (works for all sclaes)
    label = \(datetime) ymd(datetime),expand = c(0.01,0.01)) +
  #scale_x_datetime(date_breaks = "2 days", date_labels = "%Y/%m/%d",expand = c(0,0)) +
  theme(legend.text=element_text(size=17),
        text = element_text(size = 17), axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(legend.title=element_blank())

# Optional : save figure
ggsave(best.0fb0.mod.full.plot, filename='best_0fb0_full_dat.png',units='cm',width = 36, height=15)
####################################################################################################################


####################################################################################################################
# Combine figures from all four foxes into single plot ####
####################################################################################################################

myTheme <- theme(legend.text = element_text(size = 18),legend.key.size=unit(5, 'cm'))

bestpanel <- ggarrange(best.0fae.mod.plot, best.0faf.mod.plot, best.0fb1.mod.plot,
                       best.0fb0.mod.plot, 
                       labels = c("A.", "B.", "C.","D."),
                       ncol = 1, nrow = 4, common.legend = T, legend = "top") +
  myTheme

# Optional : save figure
ggsave(bestpanel, filename='bestmodels.png',units='cm',width = 28, height=35)