Last updated: 2020-12-11

Checks: 7 0

Knit directory: R_gene_analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.6.2.9000). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20200917)

The command set.seed(20200917) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 3935c51

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 3935c51. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/

Untracked files:
    Untracked:  data/Brec_R1.txt
    Untracked:  data/Brec_R2.txt
    Untracked:  data/CR15_R1.txt
    Untracked:  data/CR15_R2.txt
    Untracked:  data/CR_14_R1.txt
    Untracked:  data/CR_14_R2.txt
    Untracked:  data/KS_R1.txt
    Untracked:  data/KS_R2.txt
    Untracked:  data/NBS_PAV.txt.gz
    Untracked:  data/NLR_PAV_GD.txt
    Untracked:  data/NLR_PAV_GM.txt
    Untracked:  data/PAVs_newick.txt
    Untracked:  data/PPR1.txt
    Untracked:  data/PPR2.txt
    Untracked:  data/SNPs_newick.txt
    Untracked:  data/bac.txt
    Untracked:  data/brown.txt
    Untracked:  data/cy3.txt
    Untracked:  data/cy5.txt
    Untracked:  data/early.txt
    Untracked:  data/flowerings.txt
    Untracked:  data/foregeye.txt
    Untracked:  data/height.txt
    Untracked:  data/late.txt
    Untracked:  data/mature.txt
    Untracked:  data/motting.txt
    Untracked:  data/mvp.kin.bin
    Untracked:  data/mvp.kin.desc
    Untracked:  data/oil.txt
    Untracked:  data/pdh.txt
    Untracked:  data/protein.txt
    Untracked:  data/rust_tan.txt
    Untracked:  data/salt.txt
    Untracked:  data/seedq.txt
    Untracked:  data/seedweight.txt
    Untracked:  data/stem_termination.txt
    Untracked:  data/sudden.txt
    Untracked:  data/virus.txt
    Untracked:  data/yield.txt

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/yield_link.Rmd) and HTML (docs/yield_link.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	3935c51	Philipp Bayer	2020-12-11	wflow_publish(c(“analysis/index.Rmd”, “analysis/yield_link.Rmd”))
html	a33cd7c	Philipp Bayer	2020-12-08	Build site.
Rmd	2eec51d	Philipp Bayer	2020-12-08	wflow_publish(“analysis/yield_link.Rmd”)
html	fa16e50	Philipp Bayer	2020-11-18	Build site.
Rmd	77b5633	Philipp Bayer	2020-11-18	wflow_publish(“analysis/yield_link.Rmd”)
html	f0548a4	Philipp Bayer	2020-11-18	Build site.
Rmd	5b31dcc	Philipp Bayer	2020-11-18	wflow_publish(“analysis/yield_link.Rmd”)
html	7fda2ed	Philipp Bayer	2020-11-18	Build site.
Rmd	0e810db	Philipp Bayer	2020-11-18	More fixedp lots
html	4ed2150	Philipp Bayer	2020-11-18	Build site.
Rmd	f0de005	Philipp Bayer	2020-11-18	PLots in one row now
html	e324223	Philipp Bayer	2020-11-18	Build site.
Rmd	b8b60b1	Philipp Bayer	2020-11-18	Add StrengeJacke plots
html	fb8bbfb	Philipp Bayer	2020-11-18	Build site.
Rmd	f0a5bdb	Philipp Bayer	2020-11-18	give me back my warnings
html	adbc058	Philipp Bayer	2020-11-18	Build site.
Rmd	dc2aad4	Philipp Bayer	2020-11-18	wflow_publish(“analysis/yield_link.Rmd”)
html	359d849	Philipp Bayer	2020-11-18	Build site.
Rmd	facd3e7	Philipp Bayer	2020-11-18	wflow_publish(“analysis/yield_link.Rmd”)
html	e84338b	Philipp Bayer	2020-11-18	Build site.
Rmd	2009819	Philipp Bayer	2020-11-18	fixed lme4 syntax
html	22a48ac	Philipp Bayer	2020-11-18	Build site.
Rmd	92e2580	Philipp Bayer	2020-11-18	fixed lme4 syntax
html	262a76f	Philipp Bayer	2020-11-05	Build site.
Rmd	69e9c29	Philipp Bayer	2020-11-05	wflow_publish(“analysis/yield_link.Rmd”)
html	5ddfe2b	Philipp Bayer	2020-11-05	Build site.
Rmd	25f0f54	Philipp Bayer	2020-11-05	wflow_publish(“analysis/yield_link.Rmd”)
html	fa5c0ff	Philipp Bayer	2020-11-04	Build site.
Rmd	2d9c3db	Philipp Bayer	2020-11-04	wflow_publish(c(“analysis/index.Rmd”, “analysis/yield_link.Rmd”))
html	f34dd48	Philipp Bayer	2020-11-02	Build site.
Rmd	be2f299	Philipp Bayer	2020-11-02	wflow_publish(“analysis/yield_link.Rmd”)
html	58f8610	Philipp Bayer	2020-11-02	Build site.
Rmd	5166687	Philipp Bayer	2020-11-02	wflow_publish(“analysis/yield_link.Rmd”)
Rmd	dae157b	Philipp Bayer	2020-09-24	Update of analysis
html	dae157b	Philipp Bayer	2020-09-24	Update of analysis

knitr::opts_chunk$set(message = FALSE) 
library(tidyverse)

-- Attaching packages ------------------------------------------------------------------------------------------------------------------- tidyverse 1.3.0 --

v ggplot2 3.3.2     v purrr   0.3.4
v tibble  3.0.2     v dplyr   1.0.0
v tidyr   1.1.0     v stringr 1.4.0
v readr   1.3.1     v forcats 0.5.0

-- Conflicts ---------------------------------------------------------------------------------------------------------------------- tidyverse_conflicts() --
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(patchwork)
library(sjPlot)

Learn more about sjPlot with 'browseVignettes("sjPlot")'.

library(ggsci)
library(dabestr)

Loading required package: magrittr


Attaching package: 'magrittr'

The following object is masked from 'package:purrr':

    set_names

The following object is masked from 'package:tidyr':

    extract

library(dabestr)
library(cowplot)


********************************************************

Note: As of version 1.0.0, cowplot does not change the

  default ggplot2 theme anymore. To recover the previous

  behavior, execute:
  theme_set(theme_cowplot())

********************************************************


Attaching package: 'cowplot'

The following objects are masked from 'package:sjPlot':

    plot_grid, save_plot

The following object is masked from 'package:patchwork':

    align_plots

library(ggsignif)
library(ggforce)
library(lme4)

Loading required package: Matrix


Attaching package: 'Matrix'

The following objects are masked from 'package:tidyr':

    expand, pack, unpack

library(directlabels)
library(lmerTest)


Attaching package: 'lmerTest'

The following object is masked from 'package:lme4':

    lmer

The following object is masked from 'package:stats':

    step

library(sjPlot)
library(dotwhisker)
library(pals)
theme_set(theme_cowplot())

Data loading

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')

nbs <- read_tsv('./data/Lee.NBS.candidates.lst', col_names = c('Name', 'Class'))
nbs

# A tibble: 486 x 2
   Name                   Class
   <chr>                  <chr>
 1 UWASoyPan00953.t1      CN   
 2 GlymaLee.13G222900.1.p CN   
 3 GlymaLee.18G227000.1.p CN   
 4 GlymaLee.18G080600.1.p CN   
 5 GlymaLee.20G036200.1.p CN   
 6 UWASoyPan01876.t1      CN   
 7 UWASoyPan04211.t1      CN   
 8 GlymaLee.19G105400.1.p CN   
 9 GlymaLee.18G085100.1.p CN   
10 GlymaLee.11G142600.1.p CN   
# ... with 476 more rows

# have to remove the .t1s 
nbs$Name <- gsub('.t1','', nbs$Name)
nbs_pav_table <- pav_table %>% filter(Individual %in% nbs$Name)

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))

Warning: The `nrow` argument of `new_tibble()` can't be missing as of tibble 2.0.0.
`x` must be a scalar integer.
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

# 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))

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

# A tibble: 1,069 x 3
   `Data-storage-ID` `PI-ID`   Group   
   <chr>             <chr>     <fct>   
 1 SRR1533284        PI416890  Landrace
 2 SRR1533282        PI323576  Landrace
 3 SRR1533292        PI157421  Landrace
 4 SRR1533216        PI594615  Landrace
 5 SRR1533239        PI603336  Landrace
 6 USB-108           PI165675  Landrace
 7 HNEX-13           PI253665D Landrace
 8 USB-382           PI603549  Landrace
 9 SRR1533236        PI587552  Landrace
10 SRR1533332        PI567293  Landrace
# ... with 1,059 more rows

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

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

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

Protein

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

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

summary(lm(Protein ~ presences, data = protein_join))


Call:
lm(formula = Protein ~ presences, data = protein_join)

Residuals:
     Min       1Q   Median       3Q      Max 
-11.8479  -2.1274  -0.3336   1.9959  10.0949 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -7.98158    7.24125  -1.102    0.271    
presences    0.11786    0.01624   7.258 8.07e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.106 on 960 degrees of freedom
Multiple R-squared:  0.05203,   Adjusted R-squared:  0.05104 
F-statistic: 52.69 on 1 and 960 DF,  p-value: 8.075e-13

Seed weight

Let’s look at seed weight:

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))

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

summary(lm(wt ~ presences, data = seed_join))


Call:
lm(formula = wt ~ presences, data = seed_join)

Residuals:
     Min       1Q   Median       3Q      Max 
-12.2910  -2.8692   0.1462   2.7771  19.6962 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 91.40656   14.67990   6.227 8.28e-10 ***
presences   -0.17636    0.03298  -5.348 1.21e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.714 on 690 degrees of freedom
Multiple R-squared:  0.0398,    Adjusted R-squared:  0.0384 
F-statistic:  28.6 on 1 and 690 DF,  p-value: 1.213e-07

Oil content

And now let’s look at the oil phenotype:

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

# A tibble: 962 x 3
   names presences   Oil
   <chr>     <dbl> <dbl>
 1 AB-01       445  17.6
 2 AB-02       454  16.8
 3 BR-24       455  20.6
 4 ESS         454  20.9
 5 For         448  21  
 6 HN001       448  23.6
 7 HN002       444  18.5
 8 HN003       446  17.5
 9 HN004       442  18.9
10 HN005       440  15.5
# ... with 952 more rows

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

summary(lm(Oil ~ presences, data = oil_join))


Call:
lm(formula = Oil ~ presences, data = oil_join)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.4376  -1.9081   0.4846   2.2401   9.0361 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 118.03941    7.31646   16.13   <2e-16 ***
presences    -0.22591    0.01641  -13.77   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.139 on 960 degrees of freedom
Multiple R-squared:  0.1649,    Adjusted R-squared:  0.1641 
F-statistic: 189.6 on 1 and 960 DF,  p-value: < 2.2e-16

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

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:

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')

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:

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:

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:

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.

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:

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:

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

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

nbs_joined_groups$presences2 <- scale(nbs_joined_groups$presences, center=T, scale=T)
hist(nbs_joined_groups$presences2)

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)

basic.lm <- lm(Oil2 ~ presences2, data=oil_nbs_joined_groups)

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

Hm looks messy, you can see two groups

plot(basic.lm, which = 1)

which is confirmed by the messy line

plot(basic.lm, which = 2)

and this garbage qqplot.

So let’s build an lmer model!

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

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Oil2 ~ presences2 + (1 | Group)
   Data: oil_nbs_joined_groups

REML criterion at convergence: 1872.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.5879 -0.5672  0.0869  0.6631  3.2111 

Random effects:
 Groups   Name        Variance Std.Dev.
 Group    (Intercept) 1.3349   1.1554  
 Residual             0.4075   0.6384  
Number of obs: 951, groups:  Group, 4

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)  
(Intercept)  -0.04360    0.57867   2.99844  -0.075   0.9447  
presences2   -0.05350    0.02394 947.27006  -2.234   0.0257 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
           (Intr)
presences2 -0.004

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!

plot(mixed.lmer)

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

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')

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

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Yield2 ~ presences2 + (1 | Group)
   Data: yield_nbs_joined_groups

REML criterion at convergence: 2060.4

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1643 -0.6819  0.0316  0.6948  2.8002 

Random effects:
 Groups   Name        Variance Std.Dev.
 Group    (Intercept) 0.6466   0.8041  
 Residual             0.8600   0.9274  
Number of obs: 761, groups:  Group, 3

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)   0.23641    0.46910   1.98335   0.504 0.664692    
presences2   -0.15364    0.04172 757.46580  -3.683 0.000247 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
           (Intr)
presences2 0.025

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

plot(mixed.lmer)

qqnorm(resid(mixed.lmer))
qqline(resid(mixed.lmer))

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

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

# 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

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:

table(yield_country_nbs_joined_groups$Group)


           Wild        Landrace    Old cultivar Modern cultivar 
              0             656              33              52

And a summary histogram:

yield_country_nbs_joined_groups %>% ggplot(aes(x=presences.x, fill=Group)) + 
  geom_histogram(bins=25) +
  xlab('Yield') +
  ylab('Count') +
  facet_wrap(~Group) +
  scale_fill_manual(values = col_list) +
  theme(legend.position = "none")

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

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Yield2 ~ presences2 + (1 | Group) + (1 | Country)
   Data: yield_country_nbs_joined_groups

REML criterion at convergence: 1957

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.09429 -0.56737  0.03072  0.65680  2.89981 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.3807   0.6170  
 Group    (Intercept) 0.4178   0.6464  
 Residual             0.7614   0.8726  
Number of obs: 741, groups:  Country, 40; Group, 3

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)   
(Intercept)   0.07150    0.40194   2.28533   0.178  0.87336   
presences2   -0.11258    0.04116 726.98206  -2.735  0.00639 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
           (Intr)
presences2 0.051

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

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:

plot(mixed.lmer)

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

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

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Yield2 ~ presences2 + Group + (1 | Country)
   Data: yield_country_nbs_joined_groups

REML criterion at convergence: 1951.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1864 -0.5700  0.0305  0.6525  2.8982 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.3776   0.6145  
 Residual             0.7616   0.8727  
Number of obs: 741, groups:  Country, 40

Fixed effects:
                      Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)           -0.15265    0.13488  33.06645  -1.132  0.26588    
presences2            -0.11149    0.04117 726.45969  -2.708  0.00693 ** 
GroupOld cultivar     -0.29578    0.16164 734.95182  -1.830  0.06768 .  
GroupModern cultivar   1.00458    0.22335 360.40725   4.498 9.28e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) prsnc2 GrpOlc
presences2   0.119              
GrpOldcltvr -0.112 -0.008       
GrpMdrncltv -0.194  0.065  0.156

Non-normalised yield

Let’s see whether the ‘raw’ values perform the same.

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

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Yield ~ presences.x + (1 | Group) + (1 | Country)
   Data: yield_country_nbs_joined_groups

REML criterion at convergence: 1679.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.09429 -0.56737  0.03072  0.65680  2.89981 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.2602   0.5101  
 Group    (Intercept) 0.2856   0.5345  
 Residual             0.5205   0.7215  
Number of obs: 741, groups:  Country, 40; Group, 3

Fixed effects:
              Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)   9.011013   2.481360 677.994843   3.631 0.000303 ***
presences.x  -0.015192   0.005555 726.982171  -2.735 0.006389 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr)
presences.x -0.991

Oh, lower p-values for the intercept

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('Yield') +
      scale_color_manual(values=as.vector(isol(40)))

plot(mixed.lmer)

qqnorm(resid(mixed.lmer))
qqline(resid(mixed.lmer))

Making the breeding group fixed

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

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Yield ~ presences.x + Group + (1 | Country)
   Data: yield_country_nbs_joined_groups

REML criterion at convergence: 1675.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1864 -0.5700  0.0305  0.6525  2.8982 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.2581   0.5081  
 Residual             0.5206   0.7216  
Number of obs: 741, groups:  Country, 40

Fixed effects:
                       Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)            8.760056   2.465570 726.660933   3.553 0.000406 ***
presences.x           -0.015045   0.005556 726.459065  -2.708 0.006929 ** 
GroupOld cultivar     -0.244554   0.133648 734.951816  -1.830 0.067680 .  
GroupModern cultivar   0.830604   0.184672 360.407255   4.498 9.28e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) prsnc. GrpOlc
presences.x -0.999              
GrpOldcltvr  0.003 -0.008       
GrpMdrncltv -0.074  0.065  0.156

Oh, lower p-values for the intercept

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('Yield') +
      scale_color_manual(values=as.vector(isol(40)))

plot(mixed.lmer)

qqnorm(resid(mixed.lmer))
qqline(resid(mixed.lmer))

plot(resid(mixed.lmer))

These are the final numbers for the paper.

Plotting effect of each covariate

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

#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))

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


==================================================
                          Dependent variable:     
                     -----------------------------
                                Yield2            
--------------------------------------------------
presences2                     -0.111**           
                                (0.041)           
                                                  
GroupOld cultivar               -0.296            
                                (0.162)           
                                                  
GroupModern cultivar           1.005***           
                                (0.223)           
                                                  
Constant                        -0.153            
                                (0.135)           
                                                  
--------------------------------------------------
Observations                      741             
Log Likelihood                 -975.906           
Akaike Inf. Crit.              1963.812           
Bayesian Inf. Crit.            1991.460           
==================================================
Note:                *p<0.05; **p<0.01; ***p<0.001

library(ggeffects)
ggpredict(mixed.lmer, terms = c("presences2",  'Group'), type = "re") %>% 
   plot() +
   theme_minimal()

Let’s also plot that for non-normalised data

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('Yield')

# alright back to regular programming
library(lmerTest)
mixed.lmer <- lmer(Yield2 ~ presences2 + Group + (1|Country), data=yield_country_nbs_joined_groups)

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:

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

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Yield2 ~ presences2 + (presences2 | Group) + (1 | Country)
   Data: yield_country_nbs_joined_groups

REML criterion at convergence: 1954.8

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.09789 -0.56422  0.04471  0.67067  2.88976 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Country  (Intercept) 0.3920   0.6261       
 Group    (Intercept) 0.3858   0.6211       
          presences2  0.0310   0.1761   0.30
 Residual             0.7564   0.8697       
Number of obs: 741, groups:  Country, 40; Group, 3

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)
(Intercept)  0.02964    0.38964  2.29148   0.076    0.945
presences2  -0.22515    0.12370  1.66243  -1.820    0.235

Correlation of Fixed Effects:
           (Intr)
presences2 0.269

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

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Yield2 ~ presences2 + Group + (1 + presences2 | Country)
   Data: yield_country_nbs_joined_groups

REML criterion at convergence: 1951.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1904 -0.5715  0.0314  0.6551  2.9003 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Country  (Intercept) 0.397511 0.63048      
          presences2  0.002137 0.04623  1.00
 Residual             0.761491 0.87263      
Number of obs: 741, groups:  Country, 40

Fixed effects:
                     Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)           -0.1470     0.1387  28.1820  -1.060   0.2982    
presences2            -0.1162     0.0430  30.8335  -2.702   0.0111 *  
GroupOld cultivar     -0.3010     0.1615 733.1173  -1.863   0.0628 .  
GroupModern cultivar   1.0093     0.2192 166.8461   4.605 8.15e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) prsnc2 GrpOlc
presences2   0.334              
GrpOldcltvr -0.109 -0.017       
GrpMdrncltv -0.188  0.044  0.158
convergence code: 0
boundary (singular) fit: see ?isSingular

Oh, a significant p-value, let’s plot plot that and compare with he previous plot:

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:

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

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradient

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge: degenerate Hessian with 2 negative eigenvalues

Warning: Model failed to converge with 1 negative eigenvalue: -9.3e+04

summary(mixed.lmer)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Yield ~ presences.x + Group + (1 + presences.x | Country)
   Data: yield_country_nbs_joined_groups

REML criterion at convergence: 1675.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1863 -0.5679  0.0294  0.6532  2.8982 

Random effects:
 Groups   Name        Variance  Std.Dev.  Corr 
 Country  (Intercept) 5.191e-01 0.7204930      
          presences.x 2.402e-07 0.0004901 -0.98
 Residual             5.206e-01 0.7215284      
Number of obs: 741, groups:  Country, 40

Fixed effects:
                      Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)            8.73673    2.46492 709.51652   3.544 0.000419 ***
presences.x           -0.01499    0.00555 412.58876  -2.702 0.007186 ** 
GroupOld cultivar     -0.24420    0.13365 733.16976  -1.827 0.068096 .  
GroupModern cultivar   0.83046    0.18493 218.59478   4.491 1.15e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) prsnc. GrpOlc
presences.x -0.999              
GrpOldcltvr  0.002 -0.007       
GrpMdrncltv -0.076  0.067  0.156
convergence code: 0
unable to evaluate scaled gradient
Model failed to converge: degenerate  Hessian with 2 negative eigenvalues

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('Yield') +
      scale_color_manual(values=as.vector(isol(40)))

Quite similar, mostly downwards trajectories for each country.

And now both random slopes:

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

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Yield2 ~ presences2 + (presences2 | Group) + (1 + presences2 |  
    Country)
   Data: yield_country_nbs_joined_groups

REML criterion at convergence: 1953.7

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.11214 -0.56909  0.04459  0.66469  2.91084 

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Country  (Intercept) 0.42704  0.6535       
          presences2  0.01045  0.1022   0.81
 Group    (Intercept) 0.37523  0.6126       
          presences2  0.04201  0.2050   0.18
 Residual             0.75392  0.8683       
Number of obs: 741, groups:  Country, 40; Group, 3

Fixed effects:
            Estimate Std. Error       df t value Pr(>|t|)
(Intercept)  0.03595    0.38772  2.35869   0.093    0.933
presences2  -0.23848    0.14300  1.96304  -1.668    0.240

Correlation of Fixed Effects:
           (Intr)
presences2 0.231

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!

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)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Oil ~ presences.x + Group + (1 | Country)
   Data: oil_country_nbs_joined_groups

REML criterion at convergence: 3543.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.4112 -0.5544  0.1126  0.6528  3.0813 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.7957   0.892   
 Residual             5.0326   2.243   
Number of obs: 789, groups:  Country, 41

Fixed effects:
                      Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)           30.65976    7.34188 782.49728   4.176  3.3e-05 ***
presences.x           -0.02802    0.01654 783.56323  -1.694 0.090583 .  
GroupOld cultivar      1.12205    0.36849 784.98582   3.045 0.002404 ** 
GroupModern cultivar   1.71759    0.48119  99.35821   3.569 0.000553 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) prsnc. GrpOlc
presences.x -0.999              
GrpOldcltvr -0.002 -0.003       
GrpMdrncltv -0.075  0.066  0.153

No significance here.

tab_model(mixed.lmer, p.val='kr')

	Oil
Predictors	Estimates	CI	p
(Intercept)	30.66	16.22 – 45.10	<0.001
presences.x	-0.03	-0.06 – 0.00	0.091
Group [Old cultivar]	1.12	0.40 – 1.85	0.002
Group [Modern cultivar]	1.72	0.73 – 2.70	0.001
Random Effects
σ²	5.03
τ₀₀ _Country	0.80
ICC	0.14
N _Country	41
Observations	789
Marginal R² / Conditional R²	0.053 / 0.182

oilmod <- mixed.lmer

Protein

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)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Protein ~ presences.x + Group + (1 | Country)
   Data: protein_country_nbs_joined_groups

REML criterion at convergence: 3867.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7191 -0.6983 -0.0562  0.5934  3.5741 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.6582   0.8113  
 Residual             7.6769   2.7707  
Number of obs: 789, groups:  Country, 41

Fixed effects:
                      Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)           34.21818    9.03231 784.99934   3.788 0.000163 ***
presences.x            0.02180    0.02033 784.80098   1.072 0.283955    
GroupOld cultivar     -1.54464    0.45281 784.14871  -3.411 0.000680 ***
GroupModern cultivar  -1.06476    0.55283  74.32151  -1.926 0.057927 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) prsnc. GrpOlc
presences.x -1.000              
GrpOldcltvr -0.003 -0.001       
GrpMdrncltv -0.083  0.075  0.141

No significance here.

tab_model(mixed.lmer, p.val='kr')

	Protein
Predictors	Estimates	CI	p
(Intercept)	34.22	16.44 – 52.00	<0.001
presences.x	0.02	-0.02 – 0.06	0.285
Group [Old cultivar]	-1.54	-2.44 – -0.65	0.001
Group [Modern cultivar]	-1.06	-2.22 – 0.09	0.071
Random Effects
σ²	7.68
τ₀₀ _Country	0.66
ICC	0.08
N _Country	41
Observations	789
Marginal R² / Conditional R²	0.025 / 0.102

protmod <- mixed.lmer

Seed weight

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)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: wt ~ presences.x + Group + (1 | Country)
   Data: seed_country_nbs_joined_groups

REML criterion at convergence: 3608.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8823 -0.6388  0.0149  0.6015  4.6946 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept)  1.353   1.163   
 Residual             17.156   4.142   
Number of obs: 633, groups:  Country, 38

Fixed effects:
                       Estimate Std. Error         df t value Pr(>|t|)  
(Intercept)           15.329943  15.266389 627.894885   1.004   0.3157  
presences.x           -0.004727   0.034360 624.361120  -0.138   0.8906  
GroupOld cultivar      1.554124   0.778695 579.732427   1.996   0.0464 *
GroupModern cultivar   1.754191   0.925435   6.832254   1.896   0.1009  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) prsnc. GrpOlc
presences.x -1.000              
GrpOldcltvr -0.013  0.009       
GrpMdrncltv -0.104  0.096  0.137

Again, no significance here.

tab_model(mixed.lmer, p.val='kr')

	wt
Predictors	Estimates	CI	p
(Intercept)	15.33	-14.75 – 45.41	0.317
presences.x	-0.00	-0.07 – 0.06	0.891
Group [Old cultivar]	1.55	0.02 – 3.09	0.047
Group [Modern cultivar]	1.75	-0.23 – 3.74	0.082
Random Effects
σ²	17.16
τ₀₀ _Country	1.35
ICC	0.07
N _Country	38
Observations	633
Marginal R² / Conditional R²	0.018 / 0.090

seedmod <- mixed.lmer

The final yield model

This is the final yield model for the paper

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

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Yield ~ presences.x + Group + (1 | Country)
   Data: yield_country_nbs_joined_groups

REML criterion at convergence: 1675.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.1864 -0.5700  0.0305  0.6525  2.8982 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.2581   0.5081  
 Residual             0.5206   0.7216  
Number of obs: 741, groups:  Country, 40

Fixed effects:
                       Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)            8.760056   2.465570 726.660933   3.553 0.000406 ***
presences.x           -0.015045   0.005556 726.459065  -2.708 0.006929 ** 
GroupOld cultivar     -0.244554   0.133648 734.951816  -1.830 0.067680 .  
GroupModern cultivar   0.830604   0.184672 360.407255   4.498 9.28e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) prsnc. GrpOlc
presences.x -0.999              
GrpOldcltvr  0.003 -0.008       
GrpMdrncltv -0.074  0.065  0.156

library(ggrepel)
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('Yield') +
      scale_color_manual(values=as.vector(isol(40)))# +

      #geom_label_repel(aes(label = Country), nudge_x = 0.35, size = 4)

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('Yield') +
      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))

plot(mixed.lmer)

qqnorm(resid(mixed.lmer))
qqline(resid(mixed.lmer))

plot(resid(mixed.lmer))

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

# A tibble: 741 x 9
   names Count `PI-ID` Group presences2.x[,1] presences.y Yield Yield2.x[,1]
   <chr> <dbl> <chr>   <fct>            <dbl>       <dbl> <dbl>        <dbl>
 1 AB-01   445 PI4580~ Land~           -0.119         445  1.54      -0.774 
 2 AB-02   454 PI6037~ Land~            1.35          454  1.3       -1.06  
 3 For     448 PI5486~ Mode~            0.371         448  3.34       1.40  
 4 HN001   448 PI5186~ Mode~            0.371         448  3.54       1.64  
 5 HN005   440 PI4041~ Land~           -0.935         440  2.15      -0.0366
 6 HN006   448 PI4077~ Land~            0.371         448  2.06      -0.145 
 7 HN007   442 PI4242~ Land~           -0.609         442  1.64      -0.653 
 8 HN008   439 PI4376~ Land~           -1.10          439  2.51       0.399 
 9 HN009   445 PI4950~ Land~           -0.119         445  1.85      -0.399 
10 HN010   447 PI4689~ Land~            0.207         447  2.8        0.749 
# ... with 731 more rows, and 1 more variable: Country <chr>

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

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


==================================================
                          Dependent variable:     
                     -----------------------------
                                 Yield            
--------------------------------------------------
Count                          -0.015**           
                                (0.006)           
                                                  
GroupOld cultivar               -0.245            
                                (0.134)           
                                                  
GroupModern cultivar           0.831***           
                                (0.185)           
                                                  
Constant                       8.760***           
                                (2.466)           
                                                  
--------------------------------------------------
Observations                      741             
Log Likelihood                 -837.565           
Akaike Inf. Crit.              1687.131           
Bayesian Inf. Crit.            1714.779           
==================================================
Note:                *p<0.05; **p<0.01; ***p<0.001

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

plot_model(mixed.lmer, data=yield_country_nbs_joined_groups_renamed) +
  theme_minimal_hgrid()

plot_model(mixed.lmer, type = "pred", terms = c("Count", "Group")) +
  theme_minimal_hgrid() +
  xlab('NLR count') + 
  theme(plot.title=element_blank())

tab_model(mixed.lmer, p.val='kr')

	Yield
Predictors	Estimates	CI	p
(Intercept)	8.76	3.91 – 13.61	<0.001
Count	-0.02	-0.03 – -0.00	0.007
Group [Old cultivar]	-0.24	-0.51 – 0.02	0.068
Group [Modern cultivar]	0.83	0.46 – 1.20	<0.001
Random Effects
σ²	0.52
τ₀₀ _Country	0.26
ICC	0.33
N _Country	40
Observations	741
Marginal R² / Conditional R²	0.070 / 0.378

σ 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:

tab_model(mixed.lmer, oilmod, protmod, seedmod )

	Yield			Oil			Protein			wt
Predictors	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p	Estimates	CI	p
(Intercept)	8.76	3.93 – 13.59	<0.001	30.66	16.27 – 45.05	<0.001	34.22	16.52 – 51.92	<0.001	15.33	-14.59 – 45.25	0.315
Count	-0.02	-0.03 – -0.00	0.007
Group [Old cultivar]	-0.24	-0.51 – 0.02	0.067	1.12	0.40 – 1.84	0.002	-1.54	-2.43 – -0.66	0.001	1.55	0.03 – 3.08	0.046
Group [Modern cultivar]	0.83	0.47 – 1.19	<0.001	1.72	0.77 – 2.66	<0.001	-1.06	-2.15 – 0.02	0.054	1.75	-0.06 – 3.57	0.058
presences.x				-0.03	-0.06 – 0.00	0.090	0.02	-0.02 – 0.06	0.284	-0.00	-0.07 – 0.06	0.891
Random Effects
σ²	0.52			5.03			7.68			17.16
τ₀₀	0.26 _Country			0.80 _Country			0.66 _Country			1.35 _Country
ICC	0.33			0.14			0.08			0.07
N	40 _Country			41 _Country			41 _Country			38 _Country
Observations	741			789			789			633
Marginal R² / Conditional R²	0.070 / 0.378			0.053 / 0.182			0.025 / 0.102			0.018 / 0.090

We still need the same model for ALL genes

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

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

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Yield ~ presences.x + Group + (1 | Country)
   Data: yield_country_all_joined_groups

REML criterion at convergence: 1689.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.4280 -0.5759  0.0421  0.6532  2.7847 

Random effects:
 Groups   Name        Variance Std.Dev.
 Country  (Intercept) 0.2556   0.5056  
 Residual             0.5261   0.7253  
Number of obs: 741, groups:  Country, 40

Fixed effects:
                       Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)           3.635e+00  9.289e+00  7.219e+02   0.391   0.6956    
presences.x          -3.196e-05  1.921e-04  7.221e+02  -0.166   0.8679    
GroupOld cultivar    -2.475e-01  1.343e-01  7.351e+02  -1.843   0.0658 .  
GroupModern cultivar  8.599e-01  1.864e-01  3.577e+02   4.613 5.53e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
            (Intr) prsnc. GrpOlc
presences.x -1.000              
GrpOldcltvr -0.004  0.003       
GrpMdrncltv -0.127  0.124  0.156

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

tab_model(mixed.lmer, p.val='kr')

	Yield
Predictors	Estimates	CI	p
(Intercept)	3.64	-14.61 – 21.89	0.696
presences.x	-0.00	-0.00 – 0.00	0.868
Group [Old cultivar]	-0.25	-0.51 – 0.02	0.066
Group [Modern cultivar]	0.86	0.49 – 1.23	<0.001
Random Effects
σ²	0.53
τ₀₀ _Country	0.26
ICC	0.33
N _Country	40
Observations	741
Marginal R² / Conditional R²	0.063 / 0.369

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('Yield') +
      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))

sessionInfo()

R version 3.6.3 (2020-02-29)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17134)

Matrix products: default

locale:
[1] LC_COLLATE=English_Australia.1252  LC_CTYPE=English_Australia.1252   
[3] LC_MONETARY=English_Australia.1252 LC_NUMERIC=C                      
[5] LC_TIME=English_Australia.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] lmerTest_3.1-2         ggrepel_0.8.2          ggeffects_0.16.0      
 [4] stargazer_5.2.2        pals_1.6               dotwhisker_0.5.0      
 [7] directlabels_2020.6.17 lme4_1.1-21            Matrix_1.2-18         
[10] ggforce_0.3.1          ggsignif_0.6.0         cowplot_1.0.0         
[13] dabestr_0.3.0          magrittr_1.5           ggsci_2.9             
[16] sjPlot_2.8.6           patchwork_1.0.0        forcats_0.5.0         
[19] stringr_1.4.0          dplyr_1.0.0            purrr_0.3.4           
[22] readr_1.3.1            tidyr_1.1.0            tibble_3.0.2          
[25] ggplot2_3.3.2          tidyverse_1.3.0        workflowr_1.6.2.9000  

loaded via a namespace (and not attached):
  [1] TH.data_1.0-10      minqa_1.2.4         colorspace_1.4-1   
  [4] ellipsis_0.3.1      sjlabelled_1.1.7    rprojroot_1.3-2    
  [7] estimability_1.3    ggstance_0.3.4      parameters_0.9.0   
 [10] fs_1.5.0.9000       dichromat_2.0-0     rstudioapi_0.11    
 [13] glmmTMB_1.0.2.1     hexbin_1.28.1       farver_2.0.3       
 [16] fansi_0.4.1         mvtnorm_1.1-1       lubridate_1.7.9    
 [19] xml2_1.3.2          codetools_0.2-16    splines_3.6.3      
 [22] knitr_1.29          sjmisc_2.8.5        polyclip_1.10-0    
 [25] jsonlite_1.7.1      nloptr_1.2.1        broom_0.5.6        
 [28] dbplyr_1.4.4        effectsize_0.4.0    mapproj_1.2.7      
 [31] compiler_3.6.3      httr_1.4.2          sjstats_0.18.0     
 [34] emmeans_1.4.5       backports_1.1.10    assertthat_0.2.1   
 [37] cli_2.0.2           later_1.1.0.1       tweenr_1.0.1       
 [40] htmltools_0.5.0     tools_3.6.3         coda_0.19-3        
 [43] gtable_0.3.0        glue_1.4.2          maps_3.3.0         
 [46] Rcpp_1.0.5          cellranger_1.1.0    vctrs_0.3.1        
 [49] nlme_3.1-148        insight_0.10.0      xfun_0.17          
 [52] ps_1.3.4            rvest_0.3.5         lifecycle_0.2.0    
 [55] getPass_0.2-2       MASS_7.3-51.6       zoo_1.8-8          
 [58] scales_1.1.1        hms_0.5.3           promises_1.1.1     
 [61] sandwich_2.5-1      RColorBrewer_1.1-2  TMB_1.7.16         
 [64] yaml_2.2.1          stringi_1.5.3       bayestestR_0.7.5   
 [67] boot_1.3-25         rlang_0.4.7         pkgconfig_2.0.3    
 [70] evaluate_0.14       lattice_0.20-41     labeling_0.3       
 [73] processx_3.4.4      tidyselect_1.1.0    plyr_1.8.6         
 [76] R6_2.4.1            generics_0.0.2      multcomp_1.4-13    
 [79] DBI_1.1.0           mgcv_1.8-31         pillar_1.4.4       
 [82] haven_2.3.1         whisker_0.4         withr_2.2.0        
 [85] survival_3.2-3      performance_0.5.1   modelr_0.1.8       
 [88] crayon_1.3.4        utf8_1.1.4          rmarkdown_2.3      
 [91] grid_3.6.3          readxl_1.3.1        blob_1.2.1         
 [94] callr_3.4.4         git2r_0.27.1        reprex_0.3.0       
 [97] digest_0.6.25       xtable_1.8-4        httpuv_1.5.4       
[100] numDeriv_2016.8-1.1 munsell_0.5.0       quadprog_1.5-8

Linking yield with NLR PAV

Philipp Bayer

2020-09-22

Data loading

Linking with yield

Yield

Protein

Seed weight

Oil content

Boxplots per group

Yield

Seed weight

Oil content

Mixed modeling

Oil

Yield

Making the breeding group fixed

Adding country

Yield

Making the breeding group fixed

Non-normalised yield

Making the breeding group fixed

Plotting effect of each covariate

More complex models

Oil

Protein

Seed weight

The final yield model

We still need the same model for ALL genes