library(readxl)
library(tidyverse)
library(nlme)
library(emmeans)
library(patchwork)

theme_set(theme_bw())
palette <- c("#E66100", "#5D3A9B") #colour blind-safe

set.seed(20220823)

# Leyendo archivo original
original <- read_excel('./raw_data/data_jpaper.xlsx')

# Extrayendo datos juveniles, renombrando variables
juvenil <- original %>% 
  select(RODAL:`Coarseness (µg/m)...9`) %>% 
  rename(stand = RODAL, 
         tree = ARBOL,
         disk = `RODELA (alturas)`,
         ima = `IMA JUV`,
         length = `Length by length (mm)...5`,
         width = `Width by length (um)...6`,
         shape = `Shape by length (%)...7`,
         fines = `Fines by length (%)...8`,
         coarseness = `Coarseness (µg/m)...9`) %>% 
  mutate(type = 'Peeler core')

# Extrayendo datos maduros, renombrando variables
maduro <- original %>% 
  select(RODAL:`RODELA (alturas)`, `IMA MAD`:`Coarseness (µg/m)...15`) %>% 
  rename(stand = RODAL, 
         tree = ARBOL,
         disk = `RODELA (alturas)`,
         ima = `IMA MAD`,
         length = `Length by length (mm)...11`,
         width = `Width by length (um)...12`,
         shape = `Shape by length (%)...13`,
         fines = `Fines by length (%)...14`,
         coarseness = `Coarseness (µg/m)...15`) %>% 
  mutate(type = 'Slabwood')

# Todos los datos juntos
paper <- juvenil %>% 
  bind_rows(maduro) %>% 
  mutate(stand = factor(stand),
         tree = factor(tree),
         disk = factor(disk, levels = c('base', 'd13', '5.23')),
         disk = fct_recode(disk, Bottom = 'base', `1.3` = 'd13'),
         type = factor(type)) %>% 
  relocate(type, .after = disk)

# Estadisticas descriptivas (Table 2)
paper %>% 
  group_by(stand, disk) %>% 
  summarise_if(is.numeric, list(mean, sd)) 


# Algunos graficos para ver los datos
paper %>% 
  ggplot(aes(x = disk, y = length, color = type)) + 
  geom_jitter(alpha = 0.8, width = 0.075) + scale_colour_manual(values = palette) +
  facet_wrap(~ stand) +
  labs(x = 'Disc position', y = 'Fibre length (mm)', color = 'Wood type') +
  theme(legend.position = 'bottom') -> fiber_length_data_plot

paper %>% 
  ggplot(aes(x = disk, y = width, color = type)) +
  geom_jitter(alpha = 0.8, width = 0.075) + scale_colour_manual(values = palette) +
  facet_wrap(~ stand) +
  labs(x = 'Disc position', y = 'Fibre width (um)', color = 'Wood type') +
  theme(legend.position = 'bottom') -> fiber_width_data_plot

paper %>% 
  ggplot(aes(x = disk, y = fines, color = type)) +
  geom_jitter(alpha = 0.8, width = 0.075) + scale_colour_manual(values = palette) +
  facet_wrap(~ stand) +
  labs(x = 'Disc position', y = 'Fines (%)', color = 'Wood type') +
  theme(legend.position = 'bottom') -> fines_data_plot

paper %>% 
  ggplot(aes(x = disk, y = coarseness, color = type)) +
  geom_jitter(alpha = 0.8, width = 0.075) + scale_colour_manual(values = palette) +
  facet_wrap(~ stand) +
  labs(x = 'Disc position', y = 'Coarseness (ug/m)', color = 'Wood type') +
  theme(legend.position = 'bottom') -> coarseness_data_plot

# Todos juntos
(fiber_length_data_plot | fiber_width_data_plot) /
  (fines_data_plot | coarseness_data_plot) + plot_layout(guides = "collect") & theme(legend.position = 'bottom')

ggsave("outputs/raw_data.png", width = 22, height = 25, units = "cm")

# Ejemplo de modelo basico para una variable
##### Largo de fibra
largo_model_f <- lme(length ~ stand*disk*type , random = ~ 1|tree/disk,
                     data = paper)

# Model summary
plot(largo_model_f)
summary(largo_model_f)
anova(largo_model_f, type = 'marginal')
car::Anova(largo_model_f, type = 'III')
# Quick look at means
emmip(largo_model_f, type ~ disk | stand, CIs = TRUE) +
  scale_colour_manual(values = palette) +
  labs(x = "Disc position", y = "Predicted fibre length", 
       colour = "Wood type") -> fiber_length_means_plot
emmeans(largo_model_f, pairwise ~ disk | type)

##### Ancho de fibra
ancho_model_f <- lme(width ~ stand*disk*type , random = ~ 1|tree/disk,
                     data = paper)

# Resumen de modelo: diferencias entre rodela y tipo. Interaccion rodal rodela
plot(ancho_model_f)
car::Anova(ancho_model_f, type = 'III')
summary(ancho_model_f)

# Ojeada a los promedios. 
emmip(ancho_model_f, type ~ disk | stand, CIs = TRUE) +
  scale_colour_manual(values = palette) +
  labs(x = "Disc position", y = "Predicted fibre width", 
       colour = "Wood type") -> fiber_width_means_plot

emmeans(largo_model_f, pairwise ~ disk)
emmeans(largo_model_f, pairwise ~ type)

###### Finos
finos_model_f <- lme(fines ~ stand*disk*type , random = ~ 1|tree/disk,
                     data = paper)

# Diferencia entre tipos e interaccion tipo:rodal
plot(finos_model_f)
car::Anova(finos_model_f, type = 'III')
summary(finos_model_f)

# Ojeada a los promedios. 
emmip(finos_model_f, type ~ disk | stand, CIs = TRUE) +
  scale_colour_manual(values = palette) +
  labs(x = "Disc position", y = "Predicted fines", 
       colour = "Wood type") -> fines_means_plot
emmeans(finos_model_f, pairwise ~ type | stand)


###### Coarseness
coarse_model_f <- lme(coarseness ~ stand*disk*type , random = ~ 1|tree/disk,
                     data = paper)
# Solo diferencias de tipo
plot(coarse_model_f)
car::Anova(coarse_model_f, type = 'III')
summary(coarse_model_f)

# Ojeada a los promedios. 
emmip(coarse_model_f, type ~ disk | stand, CIs = TRUE) +
  scale_colour_manual(values = palette) +
  labs(x = "Disc position", y = "Predicted coarseness", 
       colour = "Wood type") -> coarse_means_plot

emmeans(coarse_model_f, pairwise ~ type)

# Putting all effect plots together
(fiber_length_means_plot | fiber_width_means_plot) /
  (fines_means_plot | coarse_means_plot) + plot_layout(guides = "collect") & theme(legend.position = 'bottom')

ggsave("outputs/predicted_effects.png", width = 22, height = 25, units = "cm")


#### Estimates for each level
predict(ref_grid(largo_model_f), interval = "confidence")
predict(ref_grid(ancho_model_f), interval = "confidence")
predict(ref_grid(finos_model_f), interval = "confidence")
predict(ref_grid(coarse_model_f), interval = "confidence")