longitudinal.qmd

---
title: "Longitudinal Analysis"  
description: "Mixed effect models for individual, time-dependent observations"
format:
  html: 
    df-print: kable
    code-fold: show
    code-summary: "Hide code"
    code-overflow: wrap
    toc-title: Page Contents
    toc: true
    toc-depth: 2
    toc-location: right
    number-sections: false
    html-math-method: katex
    smooth-scroll: true
editor: source
editor_options: 
  chunk_output_type: console
---
```{=html}
<style type="text/css">

body, td {
   font-size: 13pt;
}
code.r{
  font-size: 9pt;
}
pre {
  font-size: 11pt
}
</style>
```

```{r, warning=FALSE, message=FALSE, echo=FALSE}
library(tidyverse)
options(dplyr.summarise.inform = FALSE)

library(INLA)
library(here)
library(pals)

## Custom Functions
source(here("R/utilities.R"))
source_dir(here("R"))

## Read Data  
antem_df <- read_csv(here("local/bov_antemortem_2024.csv"))

# minimum date
min_date <- min(antem_df$date)

# add variables 
antem_df <- antem_df %>%
  mutate(hpdi = as.numeric(difftime(date, min_date, units = "hours")), # hours post donor inoculation 
         hpe = dpe*24, # convert dpe to hpe
         exp_type = if_else(group == "donor", "inoc", "cont"), # exposed by inoculation or direct contact
         fever = if_else(temp >= 104, "fever", "no_fever"), # temp >= 104 constitutes fever
         censor_status = if_else(group == "Group 1" | animal == "BR23-24", 0, 1), # No obs symptoms (0) in these 
  )

```

## Organize Data  
```{r}
longitud_df <- antem_df %>% 
  filter(group != "Group 1") %>%
  mutate(y_shed = if_else(nasal == 0, NA, nasal),
         y_shed = if_else(nasal == 45, 0, y_shed),
         y_serum = if_else(serum == 0, NA, serum),
         y_serum = if_else(serum == 45, 0, y_serum),
         temp_sc = as.numeric(scale(temp, scale = TRUE, center = TRUE)),
         time_step = dpe + 1)

#plot(density(longitud_df$y_serum, na.rm=T)) # ~normal
```
    
## Shedding Rates  
Run model for shedding (nasal swabs)
```{r}
return_quants <- c(0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975)

pc_prec_rw <- list(theta = list(prior="pc.prec", 
                                 param=c(1, 0.001)))

pc_prec_iid <- list(theta = list(prior="pc.prec", 
                                 param=c(3, 0.001)))

shed_form <- y_shed ~ 1 + f(time_step,
                             model = "rw1",
                             constr = TRUE,
                             scale.model = TRUE,
                             hyper=pc_prec_rw) +
                           f(animal,
                             model = "iid",
                             constr=TRUE,
                             hyper=pc_prec_iid) + 
                           f(exp_type, 
                             model = "iid",
                             constr=TRUE,
                             hyper=pc_prec_iid) +
                             temp 

shed_model <- inla(shed_form, 
                   data = longitud_df, 
                   family = "gaussian",
                   quantiles = return_quants,
                   verbose=FALSE,
                   control.compute=list(dic = TRUE, cpo = FALSE, waic = TRUE))

# summary(shed_model)
```

Plot Time Effect  
```{r fig.width=8, fig.height=6}
time_df <- shed_model$summary.random$time_step[, c(1, 4:10)]
names(time_df) <- c("DPE", "Q_0.025", "Q_0.05", "Q_0.25", "Q_0.5", "Q_0.75", "Q_0.95", "Q_0.975")

time_df$DPE <- time_df$DPE - 1

plot_time_rw(time_df)
```


Plot animals  
```{r fig.width=8, fig.height=5}
animal_df <- shed_model$summary.random$animal[, c(1, 4:10)]
names(animal_df) <- c("Animal", "Q_0.025", "Q_0.05", "Q_0.25", "Q_0.5", "Q_0.75", "Q_0.95", "Q_0.975")

plot_animals_linerange(animal_df)
```

Plot exposure type
```{r fig.width=8, fig.height=5}
exptype_df <- shed_model$summary.random$exp_type[, c(1, 4:10)]
names(exptype_df) <- c("type", "Q_0.025", "Q_0.05", "Q_0.25", "Q_0.5", "Q_0.75", "Q_0.95", "Q_0.975") 
exptype_df$type <- as.factor(exptype_df$type)

levels(exptype_df$type) <- c("Contact", "Inoculation")

plot_exposure_linerange(exptype_df)
```

Fixed effects coefficients
```{r fig.width=6, fig.height=5}
temp_post <- shed_model$marginals.fixed 

names(temp_post)
```
  
Temperature  
```{r}
plot_fixed_posterior(temp_post, model = shed_model, xlabel = "Temperature", fixed_eff = "temp")
```
  
Intercept
```{r}
plot_fixed_posterior(temp_post, model = shed_model, xlabel = "Intercept", fixed_eff = "(Intercept)")
```


## Serum  
Run model for serum
```{r}
pc_prec_rw <- list(theta = list(prior="pc.prec", 
                                 param=c(1, 0.001)))

pc_prec_iid <- list(theta = list(prior="pc.prec", 
                                 param=c(3, 0.001)))

serum_form <- y_serum ~ 1 + f(time_step,
                               model = "rw1",
                               constr = TRUE,
                               scale.model = TRUE,
                               hyper=pc_prec_rw) +
                             f(animal,
                               model = "iid",
                               constr=TRUE,
                               hyper=pc_prec_iid) + 
                             f(exp_type, 
                               model = "iid",
                               constr=TRUE,
                               hyper=pc_prec_iid) +
                               temp_sc 

serum_model <- inla(serum_form, 
                   data = longitud_df, 
                   family = "gaussian",
                   quantiles = return_quants,
                   verbose=FALSE,
                   control.compute=list(dic = TRUE, cpo = FALSE, waic = TRUE))

# summary(serum_model)
```
  
Plot Time Effect  
```{r fig.width=8, fig.height=6}
time_df <- serum_model$summary.random$time_step[, c(1, 4:10)]
names(time_df) <- c("DPE", "Q_0.025", "Q_0.05", "Q_0.25", "Q_0.5", "Q_0.75", "Q_0.95", "Q_0.975")

time_df$DPE <- time_df$DPE - 1

plot_time_rw(time_df)
```

  
Plot animals  
```{r fig.width=8, fig.height=5}
animal_df <- serum_model$summary.random$animal[, c(1, 4:10)]
names(animal_df) <- c("Animal", "Q_0.025", "Q_0.05", "Q_0.25", "Q_0.5", "Q_0.75", "Q_0.95", "Q_0.975")

plot_animals_linerange(animal_df)
```

Plot exposure type
```{r fig.width=8, fig.height=5}
exptype_df <- serum_model$summary.random$exp_type[, c(1, 4:10)]
names(exptype_df) <- c("type", "Q_0.025", "Q_0.05", "Q_0.25", "Q_0.5", "Q_0.75", "Q_0.95", "Q_0.975") 
exptype_df$type <- as.factor(exptype_df$type)

levels(exptype_df$type) <- c("Contact", "Inoculation")

plot_exposure_linerange(exptype_df)
```

Fixed effects coefficients
```{r fig.width=6, fig.height=5}
temp_post <- serum_model$marginals.fixed 

names(temp_post)
```
  
Temperature  
```{r}
plot_fixed_posterior(temp_post, model = serum_model, xlabel = "Temperature", fixed_eff = "temp_sc")
```
  
Intercept
```{r}
plot_fixed_posterior(temp_post, model = serum_model, xlabel = "Intercept", fixed_eff = "(Intercept)")
```


## Lesion Score   

Run model for lesion score
```{r}
pc_prec_rw <- list(theta = list(prior="pc.prec", 
                                 param=c(1, 0.001)))

pc_prec_iid <- list(theta = list(prior="pc.prec", 
                                 param=c(3, 0.001)))

score_form <- score ~ 1 + f(time_step,
                               model = "rw1",
                               constr = TRUE,
                               scale.model = TRUE,
                               hyper=pc_prec_rw) +
                             f(animal,
                               model = "iid",
                               constr=TRUE,
                               hyper=pc_prec_iid) + 
                             f(exp_type, 
                               model = "iid",
                               constr=TRUE,
                               hyper=pc_prec_iid) +
                               temp 

score_model <- inla(score_form, 
                   data = longitud_df, 
                   family = "gaussian",
                   quantiles = return_quants,
                   verbose=FALSE,
                   control.compute=list(dic = TRUE, cpo = FALSE, waic = TRUE))

# summary(score_model)
```
  
Plot Time Effect  
```{r fig.width=8, fig.height=6}
time_df <- score_model$summary.random$time_step[, c(1, 4:10)]
names(time_df) <- c("DPE", "Q_0.025", "Q_0.05", "Q_0.25", "Q_0.5", "Q_0.75", "Q_0.95", "Q_0.975")

time_df$DPE <- time_df$DPE - 1

plot_time_rw(time_df)
```
    
Plot animals  
```{r fig.width=8, fig.height=5}
animal_df <- score_model$summary.random$animal[, c(1, 4:10)]
names(animal_df) <- c("Animal", "Q_0.025", "Q_0.05", "Q_0.25", "Q_0.5", "Q_0.75", "Q_0.95", "Q_0.975")

plot_animals_linerange(animal_df)
```

Plot exposure type
```{r fig.width=8, fig.height=5}
exptype_df <- score_model$summary.random$exp_type[, c(1, 4:10)]
names(exptype_df) <- c("type", "Q_0.025", "Q_0.05", "Q_0.25", "Q_0.5", "Q_0.75", "Q_0.95", "Q_0.975") 
exptype_df$type <- as.factor(exptype_df$type)

levels(exptype_df$type) <- c("Contact", "Inoculation")

plot_exposure_linerange(exptype_df)
```

Fixed effects coefficients
```{r fig.width=6, fig.height=5}
temp_post <- score_model$marginals.fixed 

names(temp_post)
```
  
Temperature  
```{r}
plot_fixed_posterior(temp_post, model = score_model, xlabel = "Temperature", fixed_eff = "temp")
```
  
Intercept
```{r}
plot_fixed_posterior(temp_post, model = score_model, xlabel = "Intercept", fixed_eff = "(Intercept)")
```


## Temperature    
Run model for body temperature
```{r}
longitud_df$y_temper <- as.numeric(scale(longitud_df$temp, scale=TRUE, center=TRUE))

pc_prec_rw <- list(theta = list(prior="pc.prec", 
                                 param=c(1, 0.001)))

pc_prec_iid <- list(theta = list(prior="pc.prec", 
                                 param=c(3, 0.001)))

temper_form <- y_temper ~ -1 + f(time_step,
                                 model = "rw1",
                                 constr = TRUE,
                                 scale.model = TRUE,
                                 hyper=pc_prec_rw) +
                                f(animal,
                                 model = "iid",
                                 constr=TRUE,
                                 hyper=pc_prec_iid) + 
                                f(exp_type, 
                                 model = "iid",
                                 constr=TRUE,
                                 hyper=pc_prec_iid) +
                              score

temper_model <- inla(temper_form, 
                   data = longitud_df, 
                   family = "gaussian",
                   quantiles = return_quants,
                   verbose=FALSE,
                   control.compute=list(dic = TRUE, cpo = FALSE, waic = TRUE))

# summary(temper_model)
```
  
Plot Time Effect  
```{r fig.width=8, fig.height=6}
time_df <- temper_model$summary.random$time_step[, c(1, 4:10)]
names(time_df) <- c("DPE", "Q_0.025", "Q_0.05", "Q_0.25", "Q_0.5", "Q_0.75", "Q_0.95", "Q_0.975")

time_df$DPE <- time_df$DPE - 1

plot_time_rw(time_df)
```
      
    
Plot animals  
```{r fig.width=8, fig.height=5}
animal_df <- temper_model$summary.random$animal[, c(1, 4:10)]
names(animal_df) <- c("Animal", "Q_0.025", "Q_0.05", "Q_0.25", "Q_0.5", "Q_0.75", "Q_0.95", "Q_0.975")

plot_animals_linerange(animal_df)
```

Plot exposure type
```{r fig.width=8, fig.height=5}
exptype_df <- temper_model$summary.random$exp_type[, c(1, 4:10)]
names(exptype_df) <- c("type", "Q_0.025", "Q_0.05", "Q_0.25", "Q_0.5", "Q_0.75", "Q_0.95", "Q_0.975") 
exptype_df$type <- as.factor(exptype_df$type)

levels(exptype_df$type) <- c("Contact", "Inoculation")

plot_exposure_linerange(exptype_df)
```

Fixed effects coefficients
```{r fig.width=6, fig.height=5}
temp_post <- temper_model$marginals.fixed 

names(temp_post)
```
  
Score 
```{r}
plot_fixed_posterior(temp_post, model = temper_model, xlabel = "Score", fixed_eff = "score")
```