Schiedea sarmentosa and S. lydgatei pollen tube growth
John Powers
9/22/2020
Purpose
Test for post-pollination, pre-zygotic reproductive barriers between two plant species, Schiedea sarmentosa and S. lydgatei by measuring pollen tube growth from parental and hybrid crosses at two time points.
Response variables:
- Counts - assume a Poisson distribution for count data
- grains - the number of pollen grains on the stigma
- mid - the number of pollen tubes counted at the middle of the style
- bottom - the number of pollen tubes counted at the bottom of the style
- Proportions - assume a binomial distribution for binomial data
- ratio - the ratio of pollen tubes at the bottom / middle of the style, capped at 1 because having more pollen grains at the bottom than the middle of the style is impossible (must result from counting error). For the binomial model, the input is the number of sucesses (pollen tubes at bottom) and failures (min(bottom/mid) - bottom).
Fixed effects:
- cross - the maternal plant x the paternal plant
- hr - the time elapsed since pollination: 4 hr or 24 hr
Random effects:
- momid - maternal plant, specified by its cross number and plant number
- dadid - paternal plant, specified by its cross number and plant number
Explore data
library(knitr)
knitr::opts_chunk$set(comment=" ", warning=F, cache=T)
library(tidyverse)
library(ggsignif)
library(ggforce)
library(broom)
library(multcomp) #Multiple comparisons package - for general linear hypothesis testing
library(emmeans) #Estimated marginal means aka least square means. use old version for glmmADMB: install_version("emmeans", version = "1.3.1", repos = "http://cran.us.r-project.org")
library(glmmTMB) #generalized linear mixed models with Template Model Builder
#overwrite function to make emmeans and multcomp play nice
modelparm.glmmTMB <- function (model, coef. = function(x) fixef(x)[[component]],
vcov. = function(x) vcov(x)[[component]],
df = NULL, component="cond", ...) {
multcomp:::modelparm.default(model, coef. = coef., vcov. = vcov.,
df = df, ...)
}
signif.num <- function(x) {
symnum(x, corr = FALSE, na = FALSE, legend = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
symbols = c("***", "**", "*", ".", " "))
}Read data
#setwd("~/MyDocs/MEGA/UCI/Schiedea/Analysis/HybridSeeds")
pops.sarm <- c("896")
pops.lydg <- c("6682")
tubes <- read_csv("pollentubes_lydsar.csv") %>%
filter(keep == 1, hr != 2) %>% #exclude the resamples of pollinations that gave no tubes, and counts at 2 hrs
mutate(
hr = factor(paste0(hr,"hr"), levels=c("4hr","24hr")),
ratioraw = bottom/mid,
adjmid = ifelse(bottom > mid, bottom, mid),
ratio = bottom/adjmid,
justmid = mid - bottom,
momid = paste(mompop, momplant),
dadid = paste(dadpop, dadplant),
momcross = ifelse(mompop %in% pops.sarm, "S", "L"),
dadcross = ifelse(dadpop %in% pops.sarm, "S", "L"),
cross = factor(paste0(momcross,"x",dadcross), levels=c("LxL","LxS","SxS","SxL")),
crosshr = factor(paste0(cross, hr), levels = paste(rep(levels(cross),each=2), levels(hr), sep=""))) %>%
mutate_if(is.character, as.factor)
str(tubes) tibble [180 × 22] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
$ line : num [1:180] 2 3 4 5 6 7 8 9 10 11 ...
$ emasculation_date: Date[1:180], format: "2020-07-10" "2020-07-10" ...
$ pollination_date : Date[1:180], format: "2020-07-13" "2020-07-13" ...
$ mompop : num [1:180] 6682 6682 6682 6682 6682 ...
$ momplant : Factor w/ 11 levels "12-3","122-2",..: 10 10 10 10 10 10 10 10 10 10 ...
$ dadpop : num [1:180] 6682 6682 6682 6682 6682 ...
$ dadplant : Factor w/ 15 levels "0107-1-6","12-2",..: 13 13 13 13 12 12 12 12 5 5 ...
$ hr : Factor w/ 2 levels "4hr","24hr": 1 1 2 2 1 1 2 2 1 1 ...
$ vial : num [1:180] 1 1 2 2 3 3 4 4 5 5 ...
$ mid : num [1:180] 12 4 13 3 5 7 12 7 9 9 ...
$ bottom : num [1:180] 12 3 11 4 5 7 12 7 6 5 ...
$ keep : num [1:180] 1 1 1 1 1 1 1 1 1 1 ...
$ ratioraw : num [1:180] 1 0.75 0.846 1.333 1 ...
$ adjmid : num [1:180] 12 4 13 4 5 7 12 7 9 9 ...
$ ratio : num [1:180] 1 0.75 0.846 1 1 ...
$ justmid : num [1:180] 0 1 2 -1 0 0 0 0 3 4 ...
$ momid : Factor w/ 11 levels "6682 204-1","6682 212-1",..: 7 7 7 7 7 7 7 7 7 7 ...
$ dadid : Factor w/ 15 levels "6682 204-1","6682 209-1",..: 8 8 8 8 7 7 7 7 14 14 ...
$ momcross : Factor w/ 2 levels "L","S": 1 1 1 1 1 1 1 1 1 1 ...
$ dadcross : Factor w/ 2 levels "L","S": 1 1 1 1 1 1 1 1 2 2 ...
$ cross : Factor w/ 4 levels "LxL","LxS","SxS",..: 1 1 1 1 1 1 1 1 2 2 ...
$ crosshr : Factor w/ 8 levels "LxL4hr","LxL24hr",..: 1 1 2 2 1 1 2 2 3 3 ...
- attr(*, "spec")=
.. cols(
.. line = col_double(),
.. emasculation_date = col_date(format = ""),
.. pollination_date = col_date(format = ""),
.. mompop = col_double(),
.. momplant = col_character(),
.. dadpop = col_double(),
.. dadplant = col_character(),
.. hr = col_double(),
.. vial = col_double(),
.. mid = col_double(),
.. bottom = col_double(),
.. keep = col_double()
.. )tube.contrasts <- c("SxS4hr - SxS24hr == 0",
"SxL4hr - SxL24hr == 0",
"SxS4hr - SxL4hr == 0",
"SxS24hr - SxL24hr == 0",
"LxS4hr - LxS24hr == 0",
"LxL4hr - LxL24hr == 0",
"LxS4hr - LxL4hr == 0",
"LxS24hr - LxL24hr == 0",
"SxS4hr - LxL4hr == 0",
"SxS24hr - LxL24hr == 0",
"SxS4hr+SxL4hr - LxL4hr-LxS4hr == 0",
"SxS24hr+SxL24hr - LxL24hr-LxS24hr == 0")Inventory of attempted crosses
with(tubes, table(cross, hr)) hr
cross 4hr 24hr
LxL 28 28
LxS 28 28
SxS 16 16
SxL 18 18with(tubes, table(cross, momid)) momid
cross 6682 204-1 6682 212-1 6682 222-1 6682 229-2 6682 233-2 6682 234-1
LxL 8 8 8 8 8 8
LxS 8 8 8 8 8 8
SxS 0 0 0 0 0 0
SxL 0 0 0 0 0 0
momid
cross 6682 235-2 896 12-3 896 122-2 896 123-16 896 8-2
LxL 8 0 0 0 0
LxS 8 0 0 0 0
SxS 0 8 8 8 8
SxL 0 12 8 8 8with(tubes, table(dadid, momid)) momid
dadid 6682 204-1 6682 212-1 6682 222-1 6682 229-2 6682 233-2
6682 204-1 0 4 4 0 0
6682 209-1 0 0 0 0 0
6682 212-1 0 0 0 0 0
6682 217-1 0 0 4 0 0
6682 222-1 0 4 0 4 4
6682 229-2 0 0 0 0 0
6682 233-2 4 0 0 4 0
6682 234-1 0 0 0 0 4
6682 235-2 4 0 0 0 0
896 0107-1-6 0 4 4 0 0
896 12-2 4 0 0 4 4
896 12-3 0 4 0 0 4
896 14-3 0 0 4 0 0
896 17-5 4 0 0 4 0
896 9 2-1 0 0 0 0 0
momid
dadid 6682 234-1 6682 235-2 896 12-3 896 122-2 896 123-16 896 8-2
6682 204-1 4 0 4 0 4 0
6682 209-1 0 0 4 0 4 0
6682 212-1 0 0 4 0 0 0
6682 217-1 0 0 0 0 0 0
6682 222-1 0 0 0 4 0 0
6682 229-2 0 0 0 0 0 4
6682 233-2 0 4 0 4 0 0
6682 234-1 0 4 0 0 0 0
6682 235-2 4 0 0 0 0 4
896 0107-1-6 0 4 4 0 4 0
896 12-2 0 0 0 0 0 0
896 12-3 0 0 0 4 0 0
896 14-3 4 0 4 0 0 0
896 17-5 4 4 0 4 4 4
896 9 2-1 0 0 0 0 0 4Plot raw data
Pollen tubes
mytheme <- theme(legend.text=element_text(size=rel(1)), legend.position="right", axis.text = element_text(colour="black", size=rel(1)), text=element_text(size=14), axis.ticks.x = element_blank())
timecolors <- scale_color_brewer(palette=3, type="qual", labels=c("4 hr","24 hr"))
crosslabs <- gsub("x", " x ",levels(tubes$cross), fixed=T)
(tubeplot.mid <-
ggplot(tubes, aes(x=cross, y=mid, color=hr)) +
labs(x="", y="Pollen tubes at middle of style", color="Time") +
geom_violin(position=position_dodge(), draw_quantiles=c(0.25,0.5,0.75), size=0.5, bw=1) +
geom_sina(bw=0.1) +
scale_x_discrete(labels=crosslabs) +
theme_classic() + mytheme + timecolors)
(tubeplot.bottom <-
ggplot(tubes, aes(x=cross, y=bottom, color=hr)) +
labs(x="", y="Pollen tubes at bottom of style", color="Time") +
geom_violin(position=position_dodge(), draw_quantiles=c(0.25,0.5,0.75), size=0.5, bw=1) +
geom_sina(bw=0.1)+
scale_x_discrete(labels=crosslabs) +
theme_classic() + mytheme + timecolors)
(tubeplot.ratio <-
ggplot(tubes, aes(x=cross, y=ratio, color=hr)) +
labs(x="", y="Ratio of pollen tubes at bottom : middle of style",color="Time") +
geom_violin(position=position_dodge(), draw_quantiles=c(0.25,0.5,0.75), size=0.5, bw=0.1) +
geom_sina(bw=0.1)+
scale_y_continuous(expand = expansion(add=c(0,0.01)), breaks = scales::pretty_breaks(n = 5)) +
scale_x_discrete(labels=crosslabs) +
theme_classic() + mytheme + timecolors)
(tubeplot.justmid <-
ggplot(tubes, aes(x=cross, y=justmid, color=hr)) +
labs(x="", y="Difference of pollen tubes at middle and bottom of style",color="Time") +
geom_violin(position=position_dodge(), draw_quantiles=c(0.25,0.5,0.75), size=0.5, bw=1) +
geom_sina(bw=1) +
scale_x_discrete(labels=crosslabs) +
theme_classic() + mytheme + timecolors)
(tubeplot.mid.bottom <-
ggplot(tubes, aes(x=bottom, y=mid)) +
labs(x="Pollen tubes at bottom of style", y="Pollen tubes at middle of style", color="Cross", linetype="Time") +
geom_abline(slope=1, size=1) +
geom_point(color="grey") +
stat_ellipse(aes(linetype=hr, color=cross), size=1) +
coord_fixed(xlim=c(0,NA),ylim=c(0,NA))+
scale_color_discrete(labels=crosslabs)+ scale_linetype_discrete(labels=c("4 hr","24 hr"))+
theme_classic()+theme(legend.text=element_text(size=rel(1)), legend.position="right", axis.text = element_text(colour="black", size=rel(1)), text=element_text(size=14)))
Models
Pollen tubes at middle of style
#Quasi-Poisson (NB1), log link, zero inflation
mod <- glmmTMB(mid~cross*hr + (1|momid) + (1|dadid), data=tubes, family="poisson")
summary(mod) Family: poisson ( log )
Formula: mid ~ cross * hr + (1 | momid) + (1 | dadid)
Data: tubes
AIC BIC logLik deviance df.resid
952.9 984.8 -466.4 932.9 170
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 0.09590 0.3097
dadid (Intercept) 0.01766 0.1329
Number of obs: 180, groups: momid, 11; dadid, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.93749 0.14710 13.171 < 2e-16 ***
crossLxS -0.11406 0.13050 -0.874 0.38212
crossSxS -0.13132 0.24242 -0.542 0.58802
crossSxL -0.63984 0.24295 -2.634 0.00845 **
hr24hr 0.19976 0.09561 2.089 0.03667 *
crossLxS:hr24hr -0.01932 0.13967 -0.138 0.88999
crossSxS:hr24hr -0.49900 0.17143 -2.911 0.00361 **
crossSxL:hr24hr -0.02207 0.18616 -0.119 0.90563
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1crosshr.x <- data.frame(crosshr=levels(tubes$crosshr), x=rep(1:4, each=2)+rep(c(-0.25, 0.25), 4))
cont <- glmmTMB(mid~crosshr + (1|momid) + (1|dadid), data=tubes, family="poisson") %>%
glht(linfct = mcp(crosshr=tube.contrasts)) %>%
summary(test=adjusted(type="fdr")) %>% tidy %>%
separate(contrast, c("lvl1", "lvl2"), sep = " - ", extra="merge") %>%
mutate(signif = adj.p.value < 0.05,
stars = signif.num(adj.p.value),
height = c(rep(c(1,1,2,3),2),rep(4,4)),
x1 = crosshr.x[match(lvl1, crosshr.x$crosshr),"x"],
x2 = crosshr.x[match(lvl2, crosshr.x$crosshr),"x"])
kable(cont[,-c(1,4, 7, 9,11:13)], digits=c(0,0,3,3,3,8), caption="Contrasts")| lvl1 | lvl2 | estimate | std.error | adj.p.value | stars |
|---|---|---|---|---|---|
| SxS4hr | SxS24hr | 0.299 | 0.142 | 0.088 | . |
| SxL4hr | SxL24hr | -0.178 | 0.160 | 0.372 | |
| SxS4hr | SxL4hr | 0.509 | 0.171 | 0.036 | * |
| SxS24hr | SxL24hr | 0.032 | 0.174 | 0.855 | |
| LxS4hr | LxS24hr | -0.180 | 0.102 | 0.153 | |
| LxL4hr | LxL24hr | -0.200 | 0.096 | 0.088 | . |
| LxS4hr | LxL4hr | -0.114 | 0.131 | 0.458 | |
| LxS24hr | LxL24hr | -0.133 | 0.123 | 0.372 | |
| SxS4hr | LxL4hr | -0.131 | 0.242 | 0.641 | |
| SxS24hr | LxL24hr | -0.630 | 0.247 | 0.042 | * |
| SxS4hr + SxL4hr | LxL4hr - LxS4hr | -0.657 | 0.435 | 0.225 | |
| SxS24hr + SxL24hr | LxL24hr - LxS24hr | -1.159 | 0.434 | 0.042 | * |
timefills <- scale_fill_brewer(palette=3, type="qual", labels=c("24 hr","4 hr"), direction=-1)
cld.mod <- ref_grid(mod) %>%
emmeans(~ cross*hr) %>%
summary %>%
mutate(response = exp(emmean), uSE = exp(emmean+SE), lSE = exp(emmean-SE))
kable(cld.mod[, -c(3:7)], digits=1, caption = "Estimated marginal means")| cross | hr | response | uSE | lSE |
|---|---|---|---|---|
| LxL | 4hr | 6.9 | 8.0 | 6.0 |
| LxS | 4hr | 6.2 | 7.2 | 5.3 |
| SxS | 4hr | 6.1 | 7.4 | 5.0 |
| SxL | 4hr | 3.7 | 4.5 | 3.0 |
| LxL | 24hr | 8.5 | 9.8 | 7.3 |
| LxS | 24hr | 7.4 | 8.6 | 6.4 |
| SxS | 24hr | 4.5 | 5.5 | 3.7 |
| SxL | 24hr | 4.4 | 5.3 | 3.6 |
(mid.emplot <- ggplot(cld.mod, aes(y=response, x=cross, fill=hr)) +
geom_signif(y_position=cont$height+10, xmin = cont$x1, xmax = cont$x2, annotations=cont$stars, vjust=0.5) +
geom_col(position=position_dodge2()) +
geom_linerange(aes(ymin=lSE, ymax=uSE), position=position_dodge2(0.9)) +
labs(x="", y="Pollen tubes at middle of style",fill="Time") +
scale_y_continuous(limits=c(0,13), expand = expansion(add=c(0,0.5)), breaks = scales::pretty_breaks(n = 8)) +
scale_x_discrete(labels=crosslabs) +
theme_classic() + mytheme + timefills)
Pollen tubes at bottom of style
#Quasi-Poisson (NB1), log link, zero inflation
mod <- glmmTMB(bottom~cross*hr + (1|momid) + (1|dadid), data=tubes, family="poisson")
summary(mod) Family: poisson ( log )
Formula: bottom ~ cross * hr + (1 | momid) + (1 | dadid)
Data: tubes
AIC BIC logLik deviance df.resid
981.2 1013.1 -480.6 961.2 170
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 0.19911 0.4462
dadid (Intercept) 0.04544 0.2132
Number of obs: 180, groups: momid, 11; dadid, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.1814 0.2141 5.519 3.42e-08 ***
crossLxS -0.5992 0.2107 -2.844 0.00446 **
crossSxS 0.4859 0.3389 1.434 0.15166
crossSxL -0.1144 0.3305 -0.346 0.72925
hr24hr 0.8704 0.1203 7.234 4.71e-13 ***
crossLxS:hr24hr 0.5206 0.1939 2.684 0.00727 **
crossSxS:hr24hr -1.2310 0.1919 -6.416 1.40e-10 ***
crossSxL:hr24hr -0.6656 0.2049 -3.249 0.00116 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cont <- glmmTMB(bottom~crosshr + (1|momid) + (1|dadid), data=tubes, family="poisson") %>%
glht(linfct = mcp(crosshr=tube.contrasts)) %>%
summary(test=adjusted(type="fdr")) %>% tidy %>%
separate(contrast, c("lvl1", "lvl2"), sep = " - ", extra="merge") %>%
mutate(signif = adj.p.value < 0.05,
stars = signif.num(adj.p.value),
height = c(rep(c(1,1,2,3),2),rep(4,4)),
x1 = crosshr.x[match(lvl1, crosshr.x$crosshr),"x"],
x2 = crosshr.x[match(lvl2, crosshr.x$crosshr),"x"])
kable(cont[,-c(1,4, 7, 9,11:13)], digits=c(0,0,3,3,3,8), caption="Contrasts")| lvl1 | lvl2 | estimate | std.error | adj.p.value | stars |
|---|---|---|---|---|---|
| SxS4hr | SxS24hr | 0.361 | 0.149 | 0.032 | * |
| SxL4hr | SxL24hr | -0.205 | 0.166 | 0.260 | |
| SxS4hr | SxL4hr | 0.600 | 0.203 | 0.012 | * |
| SxS24hr | SxL24hr | 0.035 | 0.206 | 0.865 | |
| LxS4hr | LxS24hr | -1.391 | 0.152 | 0.000 | *** |
| LxL4hr | LxL24hr | -0.870 | 0.120 | 0.000 | *** |
| LxS4hr | LxL4hr | -0.599 | 0.211 | 0.013 | * |
| LxS24hr | LxL24hr | -0.079 | 0.157 | 0.672 | |
| SxS4hr | LxL4hr | 0.486 | 0.339 | 0.202 | |
| SxS24hr | LxL24hr | -0.745 | 0.336 | 0.046 | * |
| SxS4hr + SxL4hr | LxL4hr - LxS4hr | 0.971 | 0.614 | 0.171 | |
| SxS24hr + SxL24hr | LxL24hr - LxS24hr | -1.446 | 0.598 | 0.032 | * |
cld.mod <- ref_grid(mod) %>%
emmeans(~ cross*hr) %>%
summary %>%
mutate(response = exp(emmean), uSE = exp(emmean+SE), lSE = exp(emmean-SE))
kable(cld.mod[, -c(3:7)], digits=1, caption = "Estimated marginal means")| cross | hr | response | uSE | lSE |
|---|---|---|---|---|
| LxL | 4hr | 3.3 | 4.0 | 2.6 |
| LxS | 4hr | 1.8 | 2.3 | 1.4 |
| SxS | 4hr | 5.3 | 6.9 | 4.1 |
| SxL | 4hr | 2.9 | 3.8 | 2.2 |
| LxL | 24hr | 7.8 | 9.5 | 6.4 |
| LxS | 24hr | 7.2 | 8.8 | 5.9 |
| SxS | 24hr | 3.7 | 4.8 | 2.8 |
| SxL | 24hr | 3.6 | 4.7 | 2.7 |
(bottom.emplot <- ggplot(cld.mod, aes(y=response, x=cross, fill=hr)) +
geom_signif(y_position=cont$height+10, xmin = cont$x1, xmax = cont$x2, annotations=cont$stars, vjust=0.5) +
geom_col(position=position_dodge2()) +
geom_linerange(aes(ymin=lSE, ymax=uSE), position=position_dodge2(0.9)) +
labs(x="", y="Pollen tubes at bottom of style",fill="Time") +
scale_y_continuous(limits=c(0,13), expand = expansion(add=c(0,0.5)), breaks = scales::pretty_breaks(n = 8)) +
scale_x_discrete(labels=crosslabs) +
theme_classic() + mytheme + timefills)
Pollen tubes at bottom : middle of style
#Beta binomial - for overdispersed binomial data. Input is the number of successes and failures.
mod <- glmmTMB(cbind(bottom, adjmid-bottom)~cross*hr + (1|momid) + (1|dadid), data=tubes, family="binomial")
summary(mod) Family: binomial ( logit )
Formula:
cbind(bottom, adjmid - bottom) ~ cross * hr + (1 | momid) + (1 | dadid)
Data: tubes
AIC BIC logLik deviance df.resid
499.6 531.6 -239.8 479.6 170
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 1.6174 1.2718
dadid (Intercept) 0.5909 0.7687
Number of obs: 180, groups: momid, 11; dadid, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.2770 0.5907 -0.469 0.63916
crossLxS -1.2789 0.5320 -2.404 0.01622 *
crossSxS 2.5537 0.9870 2.587 0.00967 **
crossSxL 1.0687 0.9083 1.177 0.23935
hr24hr 3.0171 0.3191 9.455 < 2e-16 ***
crossLxS:hr24hr 2.3357 0.5456 4.281 1.86e-05 ***
crossSxS:hr24hr -3.2107 0.5736 -5.598 2.17e-08 ***
crossSxL:hr24hr -1.5243 0.6147 -2.480 0.01315 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1cont <- glmmTMB(cbind(bottom, adjmid-bottom)~crosshr + (1|momid) + (1|dadid), data=tubes, family="binomial") %>%
glht(linfct = mcp(crosshr=tube.contrasts)) %>%
summary(test=adjusted(type="fdr")) %>% tidy %>%
separate(contrast, c("lvl1", "lvl2"), sep = " - ", extra="merge") %>%
mutate(signif = adj.p.value < 0.05,
stars = signif.num(adj.p.value),
height = c(rep(c(1,1,2,3),2),4,5,NA,NA),
x1 = crosshr.x[match(lvl1, crosshr.x$crosshr),"x"],
x2 = crosshr.x[match(lvl2, crosshr.x$crosshr),"x"])
kable(cont[,-c(1,4, 7, 9,11:13)], digits=c(0,0,3,3,3,8), caption="Contrasts")| lvl1 | lvl2 | estimate | std.error | adj.p.value | stars |
|---|---|---|---|---|---|
| SxS4hr | SxS24hr | 0.194 | 0.475 | 0.746 | |
| SxL4hr | SxL24hr | -1.493 | 0.527 | 0.014 | * |
| SxS4hr | SxL4hr | 1.485 | 0.652 | 0.039 | * |
| SxS24hr | SxL24hr | -0.201 | 0.691 | 0.771 | |
| LxS4hr | LxS24hr | -5.353 | 0.488 | 0.000 | *** |
| LxL4hr | LxL24hr | -3.017 | 0.319 | 0.000 | *** |
| LxS4hr | LxL4hr | -1.279 | 0.532 | 0.032 | * |
| LxS24hr | LxL24hr | 1.057 | 0.615 | 0.129 | |
| SxS4hr | LxL4hr | 2.554 | 0.987 | 0.023 | * |
| SxS24hr | LxL24hr | -0.657 | 1.013 | 0.620 | |
| SxS4hr + SxL4hr | LxL4hr - LxS4hr | 4.901 | 1.705 | 0.014 | * |
| SxS24hr + SxL24hr | LxL24hr - LxS24hr | -2.169 | 1.769 | 0.293 |
cld.mod <- ref_grid(mod) %>%
emmeans(~ cross*hr) %>%
summary %>%
mutate(response = plogis(emmean), uSE = plogis(emmean+SE), lSE = plogis(emmean-SE))
kable(cld.mod[, -c(3:7)], digits=2, caption = "Estimated marginal means")| cross | hr | response | uSE | lSE |
|---|---|---|---|---|
| LxL | 4hr | 0.43 | 0.58 | 0.30 |
| LxS | 4hr | 0.17 | 0.28 | 0.10 |
| SxS | 4hr | 0.91 | 0.96 | 0.82 |
| SxL | 4hr | 0.69 | 0.83 | 0.50 |
| LxL | 24hr | 0.94 | 0.97 | 0.89 |
| LxS | 24hr | 0.98 | 0.99 | 0.96 |
| SxS | 24hr | 0.89 | 0.95 | 0.78 |
| SxL | 24hr | 0.91 | 0.96 | 0.82 |
(ratio.emplot <- ggplot(cld.mod, aes(y=response, x=cross, fill=hr)) +
geom_signif(y_position=cont$height/15+1, xmin = cont$x1, xmax = cont$x2, annotations=cont$stars, vjust=0.5) +
geom_col(position=position_dodge2()) +
geom_linerange(aes(ymin=lSE, ymax=uSE), position=position_dodge2(0.9)) +
labs(x="", y="Ratio of pollen tubes at bottom : middle of style",fill="Time") +
scale_y_continuous(expand = expansion(add=c(0,0.05)), breaks = seq(0,1, by=0.2)) +
scale_x_discrete(labels=crosslabs) +
theme_classic() + mytheme + timefills)
Package citations
citation("glmmTMB")
To cite glmmTMB in publications use:
Mollie E. Brooks, Kasper Kristensen, Koen J. van Benthem, Arni
Magnusson, Casper W. Berg, Anders Nielsen, Hans J. Skaug, Martin
Maechler and Benjamin M. Bolker (2017). glmmTMB Balances Speed and
Flexibility Among Packages for Zero-inflated Generalized Linear Mixed
Modeling. The R Journal, 9(2), 378-400.
A BibTeX entry for LaTeX users is
@Article{,
author = {Mollie E. Brooks and Kasper Kristensen and Koen J. {van Benthem} and Arni Magnusson and Casper W. Berg and Anders Nielsen and Hans J. Skaug and Martin Maechler and Benjamin M. Bolker},
title = {{glmmTMB} Balances Speed and Flexibility Among Packages for Zero-inflated Generalized Linear Mixed Modeling},
year = {2017},
journal = {The R Journal},
url = {https://journal.r-project.org/archive/2017/RJ-2017-066/index.html},
pages = {378--400},
volume = {9},
number = {2},
}citation("emmeans")
To cite package 'emmeans' in publications use:
Russell Lenth (2020). emmeans: Estimated Marginal Means, aka
Least-Squares Means. R package version 1.5.0.
https://CRAN.R-project.org/package=emmeans
A BibTeX entry for LaTeX users is
@Manual{,
title = {emmeans: Estimated Marginal Means, aka Least-Squares Means},
author = {Russell Lenth},
year = {2020},
note = {R package version 1.5.0},
url = {https://CRAN.R-project.org/package=emmeans},
}citation("multcomp")
Please cite the multcomp package by the following reference:
Torsten Hothorn, Frank Bretz and Peter Westfall (2008). Simultaneous
Inference in General Parametric Models. Biometrical Journal 50(3),
346--363.
A BibTeX entry for LaTeX users is
@Article{,
title = {Simultaneous Inference in General Parametric Models},
author = {Torsten Hothorn and Frank Bretz and Peter Westfall},
journal = {Biometrical Journal},
year = {2008},
volume = {50},
number = {3},
pages = {346--363},
}