yield_link.Rmd

---
title: "Linking yield with NLR PAV"
author: "Philipp Bayer"
date: "2020-09-22"
output: workflowr::wflow_html
editor_options:
  chunk_output_type: console
---

```{r setup}
knitr::opts_chunk$set(message = FALSE) 
library(tidyverse)
library(patchwork)
library(sjPlot)
library(ggsci)
library(dabestr)
library(dabestr)
library(cowplot)
library(ggsignif)
library(ggforce)
library(lme4)
library(directlabels)
library(lmerTest)
library(sjPlot)
library(dotwhisker)
library(pals)
theme_set(theme_cowplot())
library(RColorBrewer)
```


# Data loading

```{r}
npg_col = pal_npg("nrc")(9)
col_list <- c(`Wild`=npg_col[8],
   Landrace = npg_col[3],
  `Old cultivar`=npg_col[2],
  `Modern cultivar`=npg_col[4])

pav_table <- read_tsv('./data/soybean_pan_pav.matrix_gene.txt.gz')
```


```{r}
nbs <- read_tsv('./data/Lee.NBS.candidates.lst', col_names = c('Name', 'Class'))
nbs
# have to remove the .t1s 
nbs$Name <- gsub('.t1','', nbs$Name)
nbs_pav_table <- pav_table %>% filter(Individual %in% nbs$Name)
```


```{r}
names <- c()
presences <- c()

for (i in seq_along(nbs_pav_table)){
  if ( i == 1) next
  thisind <- colnames(nbs_pav_table)[i]
  pavs <- nbs_pav_table[[i]]
  presents <- sum(pavs)
  names <- c(names, thisind)
  presences <- c(presences, presents)
}
nbs_res_tibb <- new_tibble(list(names = names, presences = presences))
```

```{r}
# let's make the same table for all genes too
names <- c()
presences <- c()

for (i in seq_along(pav_table)){
  if ( i == 1) next
  thisind <- colnames(pav_table)[i]
  pavs <- pav_table[[i]]
  presents <- sum(pavs)
  names <- c(names, thisind)
  presences <- c(presences, presents)
}
res_tibb <- new_tibble(list(names = names, presences = presences))
```


```{r}
groups <- read_csv('./data/Table_of_cultivar_groups.csv')
groups <- rename(groups, Group = `Group in violin table`)
groups <- groups %>% 
  mutate(Group = str_replace_all(Group, 'landrace', 'Landrace')) %>%
  mutate(Group = str_replace_all(Group, 'Old_cultivar', 'Old cultivar')) %>%
  mutate(Group = str_replace_all(Group, 'Modern_cultivar', 'Modern cultivar')) %>%
  mutate(Group = str_replace_all(Group, 'Wild-type', 'Wild'))

groups$Group <-
  factor(
    groups$Group,
    levels = c('Wild',
               'Landrace',
               'Old cultivar',
               'Modern cultivar')
  )
groups
nbs_joined_groups <-
  inner_join(nbs_res_tibb, groups, by = c('names' = 'Data-storage-ID'))
all_joined_groups <-
    inner_join(res_tibb, groups, by = c('names' = 'Data-storage-ID'))

```

# Linking with yield

Can we link the trajectory of NLR genes with the trajectory of yield across the history of soybean breeding? let's make a simple regression for now

## Yield

```{r yield_join}

yield <- read_tsv('./data/yield.txt')
yield_join <- inner_join(nbs_res_tibb, yield, by=c('names'='Line'))
```

```{r}
yield_join %>% ggplot(aes(x=presences, y=Yield)) + geom_hex() + geom_smooth() +
  xlab('NLR gene count')
```


## Protein

```{r protein_join}
protein <- read_tsv('./data/protein_phenotype.txt')
protein_join <- left_join(nbs_res_tibb, protein, by=c('names'='Line')) %>% filter(!is.na(Protein))
```

```{r}
protein_join %>% ggplot(aes(x=presences, y=Protein)) + geom_hex() + geom_smooth() +
  xlab('NLR gene count')
```


```{r}
summary(lm(Protein ~ presences, data = protein_join))
```

## Seed weight

Let's look at seed weight:

```{r seed_join}
seed_weight <- read_tsv('./data/Seed_weight_Phenotype.txt', col_names = c('names', 'wt'))
seed_join <- left_join(nbs_res_tibb, seed_weight) %>% filter(!is.na(wt))
```

```{r}
seed_join %>% filter(wt > 5) %>%  ggplot(aes(x=presences, y=wt)) + geom_hex() + geom_smooth() +
  ylab('Seed weight') +
  xlab('NLR gene count')
```


```{r}
summary(lm(wt ~ presences, data = seed_join))
```

## Oil content
And now let's look at the oil phenotype:

```{r oil_join}
oil <- read_tsv('./data/oil_phenotype.txt')
oil_join <- left_join(nbs_res_tibb, oil, by=c('names'='Line')) %>% filter(!is.na(Oil))
oil_join
```


```{r}
oil_join %>%  ggplot(aes(x=presences, y=Oil)) + geom_hex() + geom_smooth() +
  xlab('NLR gene count')
```

```{r}
summary(lm(Oil ~ presences, data = oil_join))
```

OK there are many, many outliers here. Clearly I'll have to do something fancier - for example, using the first two PCs as covariates might get rid of some of those outliers. 


# Boxplots per group

## Yield
```{r}
nbs_joined_groups %>% 
  filter(!is.na(Group)) %>% 
  inner_join(yield, by=c('names'='Line')) %>% 
  ggplot(aes(x=Group, y=Yield, fill = Group)) + 
  geom_boxplot() +
  scale_fill_manual(values = col_list) + 
  theme_minimal_hgrid() +
  theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
  geom_signif(comparisons = list(c('Old cultivar', 'Modern cultivar')), 
              map_signif_level = T) +
  guides(fill=FALSE) +
  ylab('Yield') +
  xlab('Accession group')
```

And let's check the dots:


```{r}
nbs_joined_groups %>% 
  filter(!is.na(Group)) %>% 
  inner_join(yield_join, by = 'names') %>% 
  ggplot(aes(y=presences.x, x=Yield, color=Group)) +
  geom_point() + 
  scale_color_manual(values = col_list) + 
  theme_minimal_hgrid() +
  theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +  
  ylab('NLR gene count')

```

```{r}
nbs_joined_groups %>% 
  filter(!is.na(Group)) %>% 
  inner_join(yield_join, by = 'names') %>% 
  filter(Group != 'Landrace') %>% 
  ggplot(aes(x=presences.x, y=Yield, color=Group)) +
  geom_point() + 
  scale_color_manual(values = col_list) + 
  theme_minimal_hgrid() +
  geom_smooth() +
  theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +  
  xlab('NLR gene count')

```
## Protein

protein vs. the four groups:

```{r}
nbs_joined_groups %>% 
  filter(!is.na(Group)) %>% 
  inner_join(protein, by=c('names'='Line')) %>% 
  ggplot(aes(x=Group, y=Protein, fill = Group)) + 
  geom_boxplot() +
  scale_fill_manual(values = col_list) + 
  theme_minimal_hgrid() +
  theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
  geom_signif(comparisons = list(c('Wild', 'Landrace'),
                                 c('Old cultivar', 'Modern cultivar')), 
              map_signif_level = T) +
  guides(fill=FALSE) +
  ylab('Protein') +
  xlab('Accession group')
```

## Seed weight
And seed weight:

```{r}
nbs_joined_groups %>% 
  filter(!is.na(Group)) %>% 
  inner_join(seed_join) %>% 
  ggplot(aes(x=Group, y=wt, fill = Group)) + 
  geom_boxplot() +
  scale_fill_manual(values = col_list) + 
  theme_minimal_hgrid() +
  theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
  geom_signif(comparisons = list(c('Wild', 'Landrace'),
                                 c('Old cultivar', 'Modern cultivar')), 
              map_signif_level = T) +
  guides(fill=FALSE) +
  ylab('Seed weight') +
  xlab('Accession group')
```

Wow, that's breeding!

## Oil content

And finally, Oil content:

```{r}
nbs_joined_groups %>% 
  filter(!is.na(Group)) %>% 
  inner_join(oil_join, by = 'names') %>% 
  ggplot(aes(x=Group, y=Oil, fill = Group)) + 
  geom_boxplot() +
  scale_fill_manual(values = col_list) + 
  theme_minimal_hgrid() +
  theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
  geom_signif(comparisons = list(c('Wild', 'Landrace'),
                                 c('Old cultivar', 'Modern cultivar')), 
              map_signif_level = T) +
  guides(fill=FALSE) +
  ylab('Oil content') +
  xlab('Accession group')
```

Oha, a single star. That's p < 0.05!

Let's redo the above hexplot, but also color the dots by group.

```{r}
nbs_joined_groups %>% 
  filter(!is.na(Group)) %>% 
  inner_join(oil_join, by = 'names') %>% 
  ggplot(aes(x=presences.x, y=Oil, color=Group)) +
  geom_point() + 
  scale_color_manual(values = col_list) + 
  theme_minimal_hgrid() +
  theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +  
  xlab('NLR gene count')
```

Oha, so it's the Wilds that drag this out a lot.

Let's remove them and see what it looks like:

```{r}
nbs_joined_groups %>% 
  filter(!is.na(Group)) %>% 
  inner_join(oil_join, by = 'names') %>% 
  filter(Group %in% c('Old cultivar', 'Modern cultivar')) %>% 
  ggplot(aes(x=presences.x, y=Oil, color=Group)) +
  geom_point() + 
  scale_color_manual(values = col_list) + 
  theme_minimal_hgrid() +
  theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +  
  xlab('NLR gene count') +
  geom_smooth()
```

Let's remove that one outlier:
```{r}
nbs_joined_groups %>% 
  filter(!is.na(Group)) %>% 
  inner_join(oil_join, by = 'names') %>% 
  filter(Group %in% c('Old cultivar', 'Modern cultivar')) %>% 
  filter(Oil > 13) %>% 
  ggplot(aes(x=presences.x, y=Oil, color=Group)) +
  geom_point() + 
  scale_color_manual(values = col_list) + 
  theme_minimal_hgrid() +
  theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +  
  xlab('NLR gene count') +
  geom_smooth()
```

Does the above oil content boxplot become different if we exclude the one outlier? I'd bet so


```{r}
nbs_joined_groups %>% 
  filter(!is.na(Group)) %>% 
  inner_join(oil_join, by = 'names') %>% 
  filter(names != 'USB-393') %>% 
  ggplot(aes(x=Group, y=Oil, fill = Group)) + 
  geom_boxplot() +
  scale_fill_manual(values = col_list) + 
  theme_minimal_hgrid() +
  theme(axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=12)) +
  geom_signif(comparisons = list(c('Wild', 'Landrace'),
                                 c('Old cultivar', 'Modern cultivar')), 
              map_signif_level = T) +
  guides(fill=FALSE) +
  ylab('Oil content') +
  xlab('Accession group')
```

Nope, still significantly higher in modern cultivars!

# Mixed modeling

Alright here's my hypothesis: There's a link between cultivar status (Old, Wild, Landrace, Modern), r-gene count, and yield, but it's 'hidden' by country differences. 

Great tutorial here: https://ourcodingclub.github.io/tutorials/mixed-models

So we'll have to build some lme4 models! 

## Oil


```{r}
nbs_joined_groups$presences2 <- scale(nbs_joined_groups$presences, center=T, scale=T)
hist(nbs_joined_groups$presences2)
```
```{r}
oil_nbs_joined_groups <- nbs_joined_groups %>% inner_join(oil_join, by = 'names') 
oil_nbs_joined_groups$Oil2 <- scale(oil_nbs_joined_groups$Oil, center=T, scale=T)

```

```{r}
basic.lm <- lm(Oil2 ~ presences2, data=oil_nbs_joined_groups)
```

```{r}
ggplot(oil_nbs_joined_groups, aes(x = presences2, y = Oil2)) +
  geom_point() +
  geom_smooth(method = "lm")
```

Hm looks messy, you can see two groups

```{r}
plot(basic.lm, which = 1)
```

which is confirmed by the messy line

```{r}
plot(basic.lm, which = 2)
```

and this garbage qqplot.

So let's build an lmer model!


```{r}
mixed.lmer <- lmer(Oil2 ~ presences2 + (1|Group), data=oil_nbs_joined_groups)
summary(mixed.lmer)
```

So the Variance for Group is 1.3349, that means it's 1.3349/(1.3349+0.4075) *100 = 76% of the variance is explained by the four groups!

```{r}
plot(mixed.lmer)
```

```{r}
qqnorm(resid(mixed.lmer))
qqline(resid(mixed.lmer))
```

These still look fairly bad - better than before, but the QQ plot still isn't on the line.

Let's quickly check yield too

## Yield

```{r}

yield_nbs_joined_groups <- nbs_joined_groups %>% inner_join(yield_join, by = 'names') 
yield_nbs_joined_groups$Yield2 <-scale(yield_nbs_joined_groups$Yield, center=T, scale=T)

yield_all_joined_groups <- all_joined_groups %>% inner_join(yield_join, by = 'names')
```

```{r}
mixed.lmer <- lmer(Yield2 ~ presences2 + (1|Group), data=yield_nbs_joined_groups)
summary(mixed.lmer)
```

Percentage explained by breeding group: 0.6466 / (0.6466+0.8600)*100 = 42%

```{r}
plot(mixed.lmer)
```

```{r}
qqnorm(resid(mixed.lmer))
qqline(resid(mixed.lmer))
```

:O

p-value of 0.000247 for the normalised presences while accounting for the breeding group, that's beautiful.



```{r}
ggplot(yield_nbs_joined_groups, aes(x = presences2, y = Yield2)) +
      facet_wrap(~Group, nrow=1) +   # a panel for each mountain range
      geom_point(alpha = 0.5) +
      theme_classic() +
      geom_line(data = cbind(yield_nbs_joined_groups, pred = predict(mixed.lmer)), aes(y = pred), size = 1) + 
      theme_minimal_hgrid() +
      theme(legend.position = "none") +
      xlab('Scaled and centered NLR gene count') +
      ylab('Scaled and centered yield') +
      scale_color_manual(values=as.vector(isol(40)))
```

#### Making the breeding group fixed

We have < 10 possible factors in the group, so making that fixed instead of random 
```{r}
# this doesn't work because you need at least one random effect
# mixed.lmer <- lmer(Yield2 ~ presences2 + Group, data=yield_nbs_joined_groups)
```


## Adding country

We should also add the country the plant is from as a random effect, that definitely has an influence too (perhaps a stronger one???)

### Yield

```{r}
country <- read_csv('./data/Cultivar_vs_country.csv')

names(country) <- c('names', 'PI-ID', 'Country')

yield_country_nbs_joined_groups <- yield_nbs_joined_groups %>% inner_join(country)

yield_country_all_joined_groups <- yield_all_joined_groups %>% inner_join(country)
```

I need a summary table of sample sizes:

```{r}
table(yield_country_nbs_joined_groups$Group)
```

And a summary histogram:

```{r}
yield_country_nbs_joined_groups %>% ggplot(aes(x=presences.x, fill=Group)) + 
  geom_histogram(bins=25) +
  xlab(expression(paste('Yield [Mg ', ha^-1, ']'))) +
  ylab('Count') +
  facet_wrap(~Group) +
  scale_fill_manual(values = col_list) +
  theme(legend.position = "none")

```

```{r}
mixed.lmer <- lmer(Yield2 ~ presences2 + (1|Group) + (1|Country), data=yield_country_nbs_joined_groups)
summary(mixed.lmer)
```

Nice! Yield is negatively correlated with the number of NLR genes when accounting for breeding group AND country

```{r}

ggplot(yield_country_nbs_joined_groups, aes(x = presences2, y = Yield2, colour = Country)) +
      facet_wrap(~Group, nrow=1) +   # a panel for each mountain range
      geom_point(alpha = 0.5) +
      theme_classic() +
      geom_line(data = cbind(yield_country_nbs_joined_groups, pred = predict(mixed.lmer)), aes(y = pred), size = 1) + 
      theme_minimal_hgrid() +
      theme(legend.position = "none") +
      xlab('Scaled and centered NLR gene count') +
      ylab('Scaled and centered yield') +
      scale_color_manual(values=as.vector(isol(40)))
```

Some diagnostics:
```{r}
plot(mixed.lmer)
```
```{r}
qqnorm(resid(mixed.lmer))
qqline(resid(mixed.lmer))
```

Hm, the qqplot looks slightly worse than when I use maturity group alone, interesting!

BIG DISCLAIMER: Currently, I treat country and group not as nested variables, they're independent. I think that is the way it should be in this case but I'm thinking.


#### Making the breeding group fixed

Since we have too few factors in the breeding groups we have to make that fixed, not random

```{r}
mixed.lmer <- lmer(Yield2 ~ presences2 + Group + (1|Country), data=yield_country_nbs_joined_groups)
summary(mixed.lmer)
```

#### Non-normalised yield

Let's see whether the 'raw' values perform the same.


```{r}
mixed.lmer <- lmer(Yield ~ presences.x + (1|Group) + (1|Country), data=yield_country_nbs_joined_groups)
summary(mixed.lmer)
```
Oh, lower p-values for the intercept

```{r}
ggplot(yield_country_nbs_joined_groups, aes(x = presences.x, y = Yield, colour = Country)) +
      facet_wrap(~Group, nrow=1) +   # a panel for each mountain range
      geom_point(alpha = 0.5) +
      theme_classic() +
      geom_line(data = cbind(yield_country_nbs_joined_groups, pred = predict(mixed.lmer)), aes(y = pred), size = 1) + 
      theme_minimal_hgrid() +
      theme(legend.position = "none") +
      xlab('NLR gene count') +
      xlab(expression(paste('Yield [Mg ', ha^-1, ']'))) +
      scale_color_manual(values=as.vector(isol(40)))

```

```{r}
plot(mixed.lmer)
```

```{r}
qqnorm(resid(mixed.lmer))
qqline(resid(mixed.lmer))
```

#### Making the breeding group fixed

```{r}
mixed.lmer <- lmer(Yield ~ presences.x + Group + (1|Country), data=yield_country_nbs_joined_groups)
summary(mixed.lmer)
```
Oh, lower p-values for the intercept

```{r}
ggplot(yield_country_nbs_joined_groups, aes(x = presences.x, y = Yield, colour = Country)) +
      facet_wrap(~Group, nrow=1) +   # a panel for each mountain range
      geom_point(alpha = 0.5) +
      theme_classic() +
      geom_line(data = cbind(yield_country_nbs_joined_groups, pred = predict(mixed.lmer)), aes(y = pred), size = 1) + 
      theme_minimal_hgrid() +
      theme(legend.position = "none") +
      xlab('NLR gene count') +
    xlab(expression(paste('Yield [Mg ', ha^-1, ']'))) +
      scale_color_manual(values=as.vector(isol(40)))

```

```{r}
plot(mixed.lmer)
```

```{r}
qqnorm(resid(mixed.lmer))
qqline(resid(mixed.lmer))
```

```{r}
plot(resid(mixed.lmer))
```

These are the final numbers for the paper.

#### Plotting effect of each covariate

```{r}
(re.effects <- plot_model(mixed.lmer, type = "re", show.values = TRUE))
```


```{r}
#lmerTest breaks these other packages so I better unload it and reload only lme4
detach("package:lmerTest", unload=TRUE)

yield_country_nbs_joined_groups_renamed <- yield_country_nbs_joined_groups
names(yield_country_nbs_joined_groups_renamed) <- c('names', 'presences.x', 'PI-ID', 'Group', 'presences2', 'presences.y', 'Yield', 'Yield2', 'Country')
mixed.lmer <- lmer(Yield2 ~ presences2 + Group + (1|Country), data=yield_country_nbs_joined_groups_renamed)

dwplot(mixed.lmer,
       vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2))
```

```{r}
library(stargazer)
stargazer(mixed.lmer, type = "text",
          digits = 3,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digit.separator = "")
```

```{r}
library(ggeffects)
ggpredict(mixed.lmer, terms = c("presences2",  'Group'), type = "re") %>% 
   plot() +
   theme_minimal()
```
Let's also plot that for non-normalised data

```{r}
mixed.lmer <- lmer(Yield ~ presences.x + Group + (1|Country), data=yield_country_nbs_joined_groups_renamed)
ggpredict(mixed.lmer, terms = c("presences.x",  'Group'), type = "re") %>% 
   plot() +
   theme_minimal_hgrid() +
  xlab('NLR count') +
  ylab(expression(paste('Yield [Mg ', ha^-1, ']')))
```


```{r}
# alright back to regular programming
library(lmerTest)
```


#### More complex models

If I add random slopes to either groups not much changes, I do get warnings indicating that there's not much in the data:
```{r}
mixed.lmer <- lmer(Yield2 ~ presences2 + (presences2|Group) + (1|Country), data=yield_country_nbs_joined_groups)
summary(mixed.lmer)
```


```{r}
mixed.lmer <- lmer(Yield2 ~ presences2 + Group + (1 + presences2|Country), data=yield_country_nbs_joined_groups)
summary(mixed.lmer)
```
Oh, a significant p-value, let's plot plot that and compare with he previous plot:

```{r}
ggplot(yield_country_nbs_joined_groups, aes(x = presences2, y = Yield2, colour = Country)) +
      facet_wrap(~Group, nrow=1) +   # a panel for each mountain range
      geom_point(alpha = 0.5) +
      theme_classic() +
      geom_line(data = cbind(yield_country_nbs_joined_groups, pred = predict(mixed.lmer)), aes(y = pred), size = 1) + 
      theme_minimal_hgrid() +
      theme(legend.position = "none") +
      xlab('Scaled and centered NLR gene count') +
      ylab('Scaled and centered yield') +
      scale_color_manual(values=as.vector(isol(40)))
```
Quite similar, mostly downwards trajectories for each country.

Let's do that non-normalised:
```{r}
mixed.lmer <- lmer(Yield ~ presences.x + Group + (1 + presences.x|Country), data=yield_country_nbs_joined_groups)
summary(mixed.lmer)
```


```{r}
ggplot(yield_country_nbs_joined_groups, aes(x = presences.x, y = Yield, colour = Country)) +
      facet_wrap(~Group, nrow=1) +   # a panel for each mountain range
      geom_point(alpha = 0.5) +
      theme_classic() +
      geom_line(data = cbind(yield_country_nbs_joined_groups, pred = predict(mixed.lmer)), aes(y = pred), size = 1) + 
      theme_minimal_hgrid() +
      theme(legend.position = "none") +
      xlab('NLR gene count') +
    ylab(expression(paste('Yield [Mg ', ha^-1, ']'))) +
      scale_color_manual(values=as.vector(isol(40)))
```
Quite similar, mostly downwards trajectories for each country.

And now both random slopes:

```{r}
mixed.lmer <- lmer(Yield2 ~ presences2 + (presences2|Group) + (1 + presences2|Country), data=yield_country_nbs_joined_groups)
summary(mixed.lmer)
```

```{r}
ggplot(yield_country_nbs_joined_groups, aes(x = presences2, y = Yield2, colour = Country)) +
      facet_wrap(~Group, nrow=1) +   # a panel for each mountain range
      geom_point(alpha = 0.5) +
      theme_classic() +
      geom_line(data = cbind(yield_country_nbs_joined_groups, pred = predict(mixed.lmer)), aes(y = pred), size = 1) + 
      theme_minimal_hgrid() +
      theme(legend.position = "none") +
      xlab('Scaled and centered NLR gene count') +
      ylab('Scaled and centered yield') +
      scale_color_manual(values=as.vector(isol(40)))
```
Yeah, nah

### Oil

I'm removing the wilds from the other phenotypes to make the models comparable with the yield model - the yield model uses Landrace as baseline, if I keep Wild in then the baseline is different!

```{r}
oil_country_nbs_joined_groups <- oil_nbs_joined_groups %>% inner_join(country)
oil_country_nbs_joined_groups <- oil_country_nbs_joined_groups %>% filter(Group != 'Wild')
mixed.lmer <- lmer(Oil ~ presences.x + Group + (1|Country), data=oil_country_nbs_joined_groups)
summary(mixed.lmer)
```

No significance here.

```{r}
tab_model(mixed.lmer, p.val='kr')
oilmod <- mixed.lmer
```


```{r}
table(oil_country_nbs_joined_groups$Group)

```

### Protein

```{r}
protein_nbs_joined_groups <- nbs_joined_groups %>% inner_join(protein_join, by = 'names') 
#protein_nbs_joined_groups$Protein2 <- scale(protein_nbs_joined_groups$Protein, center=T, scale=T)
protein_country_nbs_joined_groups <- protein_nbs_joined_groups %>% inner_join(country)
#protein_country_nbs_joined_groups <- rename(protein_country_nbs_joined_groups, Group=`Group in violin table`)
protein_country_nbs_joined_groups <- protein_country_nbs_joined_groups %>% filter(Group != 'Wild')

mixed.lmer <- lmer(Protein ~ presences.x + Group + (1|Country), data=protein_country_nbs_joined_groups)
summary(mixed.lmer)
```

No significance here.

```{r}
tab_model(mixed.lmer, p.val='kr')
protmod <- mixed.lmer
```

```{r}
table(protein_country_nbs_joined_groups$Group)

```

### Seed weight
```{r}
seed_nbs_joined_groups <- nbs_joined_groups %>% inner_join(seed_join, by = 'names') 
#seed_nbs_joined_groups$wt2 <- scale(seed_nbs_joined_groups$wt, center=T, scale=T)
seed_country_nbs_joined_groups <- seed_nbs_joined_groups %>% inner_join(country)
#seed_country_nbs_joined_groups <- rename(seed_country_nbs_joined_groups, Group = `Group in violin table`)
seed_country_nbs_joined_groups <- seed_country_nbs_joined_groups %>% filter(Group != 'Wild')
mixed.lmer <- lmer(wt ~ presences.x + Group + (1|Country), data=seed_country_nbs_joined_groups)
summary(mixed.lmer)
```

Again, no significance here.

```{r}
tab_model(mixed.lmer, p.val='kr')
seedmod <- mixed.lmer
```

```{r}
table(seed_country_nbs_joined_groups$Group)

```

## The final yield model

This is the final yield model for the paper

```{r}
mixed.lmer <- lmer(Yield ~ presences.x + Group + (1|Country), data=yield_country_nbs_joined_groups)
summary(mixed.lmer)
```


```{r}
ggplot(yield_country_nbs_joined_groups, aes(x = presences.x, y = Yield, colour = Country)) +
      facet_wrap(~Group, nrow=1) +   # a panel for each mountain range
      geom_point(alpha = 0.5) +
      geom_line(data = cbind(yield_country_nbs_joined_groups, pred = predict(mixed.lmer)), aes(y = pred), size = 1) + 
      theme_minimal_hgrid() +
      theme(legend.position = "none") +
      xlab('NLR gene count') +
  ylab(expression(paste('Yield [Mg ', ha^-1, ']'))) +
      scale_color_manual(values=as.vector(isol(40)))
```

```{r}

newdat <-cbind(yield_country_nbs_joined_groups, pred = predict(mixed.lmer))

newdat %>% mutate(Country2 = case_when ( Country == 'USA' ~ 'USA',
                                              Country == 'China' ~ 'China',
                                              Country == 'Korea' ~ 'Korea',
                                              Country == 'Japan' ~ 'Japan',
                                              Country == 'Russia' ~ 'Russia',
                                              TRUE ~ '')) %>% 
  ggplot(aes(x = presences.x, y = pred, colour = Country)) +
      facet_wrap(~Group, nrow=1) +   # a panel for each mountain range
      geom_line( size = 1) +
      theme_minimal_hgrid() +
      theme(legend.position = "none") +
      xlab('NLR gene count') +
  ylab(expression(paste('Yield [Mg ', ha^-1, ']'))) +
      scale_color_manual(values=as.vector(isol(40)))+
      geom_point(aes(y = Yield),alpha = 0.5) +
      geom_dl(aes(label = Country2), method='last.bumpup') +
      xlim(c(430, 480))

```

Let's just use 6 groups - 5 main countries plus the rest

```{r}
#remove that ugly yellow
mycol <- c(brewer.pal(n = 8, name = "Accent")[1:3], brewer.pal(n = 8, name = "Accent")[5:8])
```

```{r}
newdat %>% mutate(Country2 = case_when ( Country == 'USA' ~ 'USA',
                                              Country == 'China' ~ 'China',
                                              Country == 'Korea' ~ 'Korea',
                                              Country == 'Japan' ~ 'Japan',
                                              Country == 'Russia' ~ 'Russia',
                                              TRUE ~ 'Rest')) %>% 
  mutate(Country2 = factor(Country2, levels=c('China', 'Japan', 'Korea', 'Russia', 'USA', 'Rest'))) %>% 
  ggplot(aes(x = presences.x, y = pred, color = Country2)) +
      facet_wrap(~Group, nrow=1) +   # a panel for each mountain range
      geom_point(aes(y = Yield, color=Country2),alpha = 0.8, size=2) +
      geom_line(aes(y=pred, group=Country, color=Country2), size = 1.5) +
      theme_minimal_hgrid() +
      xlab('NLR gene count') +
      ylab(expression(paste('Yield [Mg ', ha^-1, ']'))) +
      scale_color_manual(values=mycol) +
      xlim(c(430, 480)) +
  labs(color = "Country")
  
```

Let's try another color scheme

```{r}
# I want only every second, stronger color of the Paired scheme
mycol <- brewer.pal(n = 12, name = "Paired")[seq(2, 12, 2)]
```

```{r}
newdat %>% mutate(Country2 = case_when ( Country == 'USA' ~ 'USA',
                                              Country == 'China' ~ 'China',
                                              Country == 'Korea' ~ 'Korea',
                                              Country == 'Japan' ~ 'Japan',
                                              Country == 'Russia' ~ 'Russia',
                                              TRUE ~ 'Rest')) %>% 
  mutate(Country2 = factor(Country2, levels=c('China', 'Japan', 'Korea', 'Russia', 'USA', 'Rest'))) %>% 
  ggplot(aes(x = presences.x, y = pred, color = Country2)) +
      facet_wrap(~Group, nrow=1) +   # a panel for each mountain range
      geom_point(aes(y = Yield, color=Country2),alpha = 0.8, size=2) +
      geom_line(aes(y=pred, group=Country, color=Country2), size = 1.5) +
      theme_minimal_hgrid() +
      xlab('NLR gene count') +
      ylab(expression(paste('Yield [Mg ', ha^-1, ']'))) +
      scale_color_manual(values=mycol) +
      xlim(c(430, 480)) +
  labs(color = "Country") +
  theme(panel.spacing = unit(0.9, "lines"),
        axis.text.x = element_text(size=10))
```

OK that's much better, nice and strong colors.

```{r}
plot(mixed.lmer)
```

```{r}
qqnorm(resid(mixed.lmer))
qqline(resid(mixed.lmer))
```

```{r}
plot(resid(mixed.lmer))
```

```{r}
detach("package:lmerTest", unload=TRUE)

yield_country_nbs_joined_groups_renamed <- yield_country_nbs_joined_groups
names(yield_country_nbs_joined_groups_renamed) <- c('names', 'Count', 'PI-ID', 'Group', 'presences2.x', 'presences.y', 'Yield', 'Yield2.x', 'Country')
mixed.lmer <- lmer(Yield ~ `Count` + Group + (1|Country), data=yield_country_nbs_joined_groups_renamed)
yield_country_nbs_joined_groups_renamed
dwplot(mixed.lmer,
       vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2))
```

```{r}
library(stargazer)
stargazer(mixed.lmer, type = "text",
          digits = 3,
          star.cutoffs = c(0.05, 0.01, 0.001),
          digit.separator = "")
```

```{r}
library(ggeffects)
ggpredict(mixed.lmer, terms = c("Count",  'Group'), type = "re") %>% 
   plot() +
   theme_minimal_hgrid() +
   xlab('NLR count') +  
   theme(plot.title=element_blank())
```

```{r}
plot_model(mixed.lmer, type = "re", sort.est = TRUE) + theme(plot.title=element_blank())

```

```{r}
plot_model(mixed.lmer, data=yield_country_nbs_joined_groups_renamed) +
  theme_minimal_hgrid() +
  theme(plot.title=element_blank())

```

```{r}
plot_model(mixed.lmer, type = "pred", terms = c("Count", "Group")) +
  theme_minimal_hgrid() +
  xlab('NLR count') + 
  ylab((expression(paste('Yield [Mg ', ha^-1, ']')))) +
  theme(plot.title=element_blank())
```

```{r}
tab_model(mixed.lmer, p.val='kr', digits=3)
```

σ measures the random effect variance I think, 0.52 is pretty good (this can easily be >1), but more useful to compare models with each other which I don't do here. 

intraclass-correlation coefficient (ICC) measures how the proportion of variance explained by the grouping structure, in this case, country


Let's compare all models in one table:

```{r}
tab_model(mixed.lmer, oilmod, protmod, seedmod, digits=3 )
```


## We still need the same model for ALL genes

What if we just see a general gene shrinkage, not just NLR-genes?

```{r}
library('lmerTest')
mixed.lmer <- lmer(Yield ~ presences.x + Group + (1|Country), data=yield_country_all_joined_groups)
summary(mixed.lmer)
```

OK good, so all genes don't have a statistically significant correlation.

```{r}
tab_model(mixed.lmer, p.val='kr', digits=3)

```

```{r}
newdat <-cbind(yield_country_all_joined_groups, pred = predict(mixed.lmer))

newdat %>% mutate(Country2 = case_when ( Country == 'USA' ~ 'USA',
                                              Country == 'China' ~ 'China',
                                              Country == 'Korea' ~ 'Korea',
                                              Country == 'Japan' ~ 'Japan',
                                              Country == 'Russia' ~ 'Russia',
                                              TRUE ~ '')) %>% 
  ggplot(aes(x = presences.x/1000, y = pred, colour = Country)) +
      facet_wrap(~Group, nrow=1) +   # a panel for each mountain range
      geom_line( size = 1) +
      theme_minimal_hgrid() +
      theme(legend.position = "none") +
      xlab('Gene count (1000s)') +
  ylab(expression(paste('Yield [Mg ', ha^-1, ']'))) +
      scale_color_manual(values=as.vector(isol(40)))+
      geom_point(aes(y = Yield),alpha = 0.5) +
      geom_dl(aes(label = Country2), method='last.bumpup')+
      xlim(c(47.900, 49.700))
```

With better color scheme:

```{r}
newdat %>% mutate(
  Country2 = case_when (
  Country == 'USA' ~ 'USA',
  Country == 'China' ~ 'China',
  Country == 'Korea' ~ 'Korea',
  Country == 'Japan' ~ 'Japan',
  Country == 'Russia' ~ 'Russia',
  TRUE ~ 'Rest'
  )
  ) %>%
  mutate(Country2 = factor(
  Country2,
  levels = c('China', 'Japan', 'Korea', 'Russia', 'USA', 'Rest')
  )) %>%
  ggplot(aes(x = presences.x / 1000, y = pred, color = Country2)) +
  facet_wrap( ~ Group, nrow = 1) +   # a panel for each mountain range
  geom_point(aes(y = Yield, color = Country2),
  alpha = 0.8,
  size = 2) +
  geom_line(aes(y = pred, group = Country, color = Country2), size = 1.5) +
  theme_minimal_hgrid() +
  xlab('Gene count') +
  ylab(expression(paste('Yield [Mg ', ha ^ -1, ']'))) +
  scale_color_manual(values = mycol) +
  xlim(c(47.900, 49.700)) +
  labs(color = "Country") +
  theme(panel.spacing = unit(0.9, "lines"),
  axis.text.x = element_text(size = 10))
```