Schiedea stellariodes and S. spergulina pollen tube growth
John Powers
2020-01-15
Purpose
Test for post-pollination, pre-zygotic reproductive barriers between two plant species, Schiedea stellariodes and S. spergulina 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.sper <- c("120289","120294","120296")
pops.stel <- c("14651","14765")
tubes <- read_csv("pollentubes_speste.csv") %>% mutate(across(ends_with("pop"), as.factor)) %>%
pivot_longer(starts_with(c("mid","bottom")), names_to=c("location","stigma"), names_sep="_", values_to="n_tubes") %>%
pivot_wider(names_from="location", values_from="n_tubes") %>% drop_na(mid) %>%
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.sper, "Sp", "St"),
dadcross = ifelse(dadpop %in% pops.sper, "Sp", "St"),
cross = factor(paste0(momcross,"x",dadcross), levels=c("StxSt","StxSp","SpxSp","SpxSt")),
crosshr = factor(paste0(cross, hr), levels = paste(rep(levels(cross),each=2), levels(hr), sep=""))) %>%
mutate_if(is.character, as.factor)
str(tubes) tibble [474 × 23] (S3: tbl_df/tbl/data.frame)
$ line : num [1:474] 1 1 1 2 2 2 3 3 3 4 ...
$ emasculation_date: Date[1:474], format: NA NA ...
$ pollination_date : Date[1:474], format: "2020-11-06" "2020-11-06" ...
$ mompop : Factor w/ 5 levels "14651","14765",..: 5 5 5 5 5 5 5 5 5 5 ...
$ momplant : Factor w/ 10 levels "11F","20J-6",..: 10 10 10 10 10 10 10 10 10 10 ...
$ dadpop : Factor w/ 5 levels "14651","14765",..: 2 2 2 2 2 2 2 2 2 2 ...
$ dadplant : Factor w/ 13 levels "12 (M)","13 (M)",..: 10 10 10 10 10 10 10 10 10 10 ...
$ hr : Factor w/ 2 levels "4hr","24hr": 1 1 1 1 1 1 2 2 2 2 ...
$ vial : num [1:474] 1 1 1 1 1 1 2 2 2 2 ...
$ count_date : Date[1:474], format: "2020-11-10" "2020-11-10" ...
$ stigma : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...
$ mid : num [1:474] 5 3 7 4 4 3 5 5 6 6 ...
$ bottom : num [1:474] 3 1 3 2 2 2 5 4 5 5 ...
$ ratioraw : num [1:474] 0.6 0.333 0.429 0.5 0.5 ...
$ adjmid : num [1:474] 5 3 7 4 4 3 5 5 6 6 ...
$ ratio : num [1:474] 0.6 0.333 0.429 0.5 0.5 ...
$ justmid : num [1:474] 2 2 4 2 2 1 0 1 1 1 ...
$ momid : Factor w/ 10 levels "120289 11F","120289 5F/M",..: 5 5 5 5 5 5 5 5 5 5 ...
$ dadid : Factor w/ 15 levels "120289 12 (M)",..: 14 14 14 14 14 14 14 14 14 14 ...
$ momcross : Factor w/ 2 levels "Sp","St": 1 1 1 1 1 1 1 1 1 1 ...
$ dadcross : Factor w/ 2 levels "Sp","St": 2 2 2 2 2 2 2 2 2 2 ...
$ cross : Factor w/ 4 levels "StxSt","StxSp",..: 4 4 4 4 4 4 4 4 4 4 ...
$ crosshr : Factor w/ 8 levels "StxSt4hr","StxSt24hr",..: 7 7 7 7 7 7 8 8 8 8 ...tube.contrasts <- c("SpxSp4hr - SpxSp24hr == 0",
"SpxSt4hr - SpxSt24hr == 0",
"SpxSp4hr - SpxSt4hr == 0",
"SpxSp24hr - SpxSt24hr == 0",
"StxSp4hr - StxSp24hr == 0",
"StxSt4hr - StxSt24hr == 0",
"StxSp4hr - StxSt4hr == 0",
"StxSp24hr - StxSt24hr == 0",
"SpxSp4hr - StxSt4hr == 0",
"SpxSp24hr - StxSt24hr == 0",
"SpxSp4hr+SpxSt4hr - StxSt4hr-StxSp4hr == 0",
"SpxSp24hr+SpxSt24hr - StxSt24hr-StxSp24hr == 0")Inventory of attempted crosses
with(tubes, table(cross, hr)) hr
cross 4hr 24hr
StxSt 78 77
StxSp 42 42
SpxSp 52 53
SpxSt 66 64with(tubes, table(cross, momid)) momid
cross 120289 11F 120289 5F/M 120294 4F 120296 8F 120296 9F 14651 4A
StxSt 0 0 0 0 0 36
StxSp 0 0 0 0 0 12
SpxSp 23 24 24 22 12 0
SpxSt 23 24 23 24 36 0
momid
cross 14765 20J-6 14765 4E 14765 6G-Cut 1 14765 9A-3
StxSt 24 35 24 36
StxSp 24 12 24 12
SpxSp 0 0 0 0
SpxSt 0 0 0 0with(tubes, table(dadid, momid)) momid
dadid 120289 11F 120289 5F/M 120294 4F 120296 8F 120296 9F 14651 4A
120289 12 (M) 0 0 12 0 12 0
120289 13 (M) 11 0 0 11 0 0
120289 7 (M) 0 12 0 0 0 12
120292 5 (M) 0 0 0 0 12 12
120292 5(M) 0 0 0 0 0 0
120294 12 (M) 0 0 0 0 0 0
120294 7 (M) 12 0 0 0 0 0
120294 8 (M) 0 12 12 11 0 0
14651 3B-Cut 1 0 12 0 0 0 12
14651 4A 0 12 0 0 0 0
14765 14B-2 0 0 0 0 0 0
14765 20J-6 11 0 0 0 12 0
14765 4E 0 0 12 0 0 12
14765 6G-Cut 1 12 0 0 12 12 0
14765 9A-3 0 0 11 12 0 0
momid
dadid 14765 20J-6 14765 4E 14765 6G-Cut 1 14765 9A-3
120289 12 (M) 12 0 0 0
120289 13 (M) 12 12 0 0
120289 7 (M) 0 0 12 0
120292 5 (M) 0 0 0 12
120292 5(M) 0 12 0 0
120294 12 (M) 0 0 12 0
120294 7 (M) 0 0 0 0
120294 8 (M) 0 0 0 12
14651 3B-Cut 1 12 0 0 12
14651 4A 0 11 0 0
14765 14B-2 12 12 0 0
14765 20J-6 0 0 12 0
14765 4E 0 0 12 0
14765 6G-Cut 1 0 0 0 12
14765 9A-3 0 0 0 0Plot 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
1866.1 1907.7 -923.0 1846.1 464
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 3.182e-02 1.784e-01
dadid (Intercept) 2.225e-10 1.492e-05
Number of obs: 474, groups: momid, 10; dadid, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.43263 0.09718 14.743 <2e-16 ***
crossStxSp -0.19820 0.10093 -1.964 0.0496 *
crossSpxSp 0.18465 0.13999 1.319 0.1872
crossSpxSt 0.10661 0.13803 0.772 0.4399
hr24hr 0.13329 0.07562 1.763 0.0780 .
crossStxSp:hr24hr 0.08847 0.13469 0.657 0.5113
crossSpxSp:hr24hr -0.11653 0.11400 -1.022 0.3067
crossSpxSt:hr24hr -0.07307 0.10959 -0.667 0.5050
---
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 |
|---|---|---|---|---|---|
| SpxSp4hr | SpxSp24hr | -0.017 | 0.085 | 0.844 | |
| SpxSt4hr | SpxSt24hr | -0.060 | 0.079 | 0.597 | |
| SpxSp4hr | SpxSt4hr | 0.078 | 0.084 | 0.597 | |
| SpxSp24hr | SpxSt24hr | 0.035 | 0.083 | 0.737 | |
| StxSp4hr | StxSp24hr | -0.222 | 0.111 | 0.234 | |
| StxSt4hr | StxSt24hr | -0.133 | 0.076 | 0.234 | |
| StxSp4hr | StxSt4hr | -0.198 | 0.101 | 0.234 | |
| StxSp24hr | StxSt24hr | -0.110 | 0.092 | 0.466 | |
| SpxSp4hr | StxSt4hr | 0.185 | 0.140 | 0.449 | |
| SpxSp24hr | StxSt24hr | 0.068 | 0.138 | 0.737 | |
| SpxSp4hr + SpxSt4hr | StxSt4hr - StxSp4hr | 0.489 | 0.261 | 0.234 | |
| SpxSp24hr + SpxSt24hr | StxSt24hr - StxSp24hr | 0.211 | 0.257 | 0.597 |
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 |
|---|---|---|---|---|
| StxSt | 4hr | 4.2 | 4.6 | 3.8 |
| StxSp | 4hr | 3.4 | 3.9 | 3.1 |
| SpxSp | 4hr | 5.0 | 5.6 | 4.6 |
| SpxSt | 4hr | 4.7 | 5.1 | 4.2 |
| StxSt | 24hr | 4.8 | 5.3 | 4.4 |
| StxSp | 24hr | 4.3 | 4.8 | 3.8 |
| SpxSp | 24hr | 5.1 | 5.7 | 4.6 |
| SpxSt | 24hr | 5.0 | 5.5 | 4.5 |
(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
1611.5 1653.1 -795.8 1591.5 464
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 0.05383 0.2320
dadid (Intercept) 0.02628 0.1621
Number of obs: 474, groups: momid, 10; dadid, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.6267 0.1447 4.332 1.48e-05 ***
crossStxSp -1.5653 0.2729 -5.736 9.68e-09 ***
crossSpxSp 0.5840 0.2064 2.829 0.004664 **
crossSpxSt 0.2094 0.1900 1.102 0.270376
hr24hr 0.7131 0.0999 7.138 9.44e-13 ***
crossStxSp:hr24hr 0.9970 0.2819 3.537 0.000405 ***
crossSpxSp:hr24hr -0.4297 0.1384 -3.106 0.001898 **
crossSpxSt:hr24hr -0.2129 0.1411 -1.509 0.131418
---
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 |
|---|---|---|---|---|---|
| SpxSp4hr | SpxSp24hr | -0.283 | 0.096 | 0.005 | ** |
| SpxSt4hr | SpxSt24hr | -0.500 | 0.100 | 0.000 | *** |
| SpxSp4hr | SpxSt4hr | 0.375 | 0.144 | 0.012 | * |
| SpxSp24hr | SpxSt24hr | 0.158 | 0.131 | 0.248 | |
| StxSp4hr | StxSp24hr | -1.710 | 0.264 | 0.000 | *** |
| StxSt4hr | StxSt24hr | -0.713 | 0.100 | 0.000 | *** |
| StxSp4hr | StxSt4hr | -1.565 | 0.273 | 0.000 | *** |
| StxSp24hr | StxSt24hr | -0.568 | 0.151 | 0.000 | *** |
| SpxSp4hr | StxSt4hr | 0.584 | 0.206 | 0.007 | ** |
| SpxSp24hr | StxSt24hr | 0.154 | 0.195 | 0.428 | |
| SpxSp4hr + SpxSt4hr | StxSt4hr - StxSp4hr | 2.359 | 0.409 | 0.000 | *** |
| SpxSp24hr + SpxSt24hr | StxSt24hr - StxSp24hr | 0.719 | 0.334 | 0.038 | * |
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 |
|---|---|---|---|---|
| StxSt | 4hr | 1.9 | 2.2 | 1.6 |
| StxSp | 4hr | 0.4 | 0.5 | 0.3 |
| SpxSp | 4hr | 3.4 | 3.9 | 2.9 |
| SpxSt | 4hr | 2.3 | 2.7 | 2.0 |
| StxSt | 24hr | 3.8 | 4.4 | 3.3 |
| StxSp | 24hr | 2.2 | 2.5 | 1.8 |
| SpxSp | 24hr | 4.5 | 5.1 | 3.9 |
| SpxSt | 24hr | 3.8 | 4.4 | 3.3 |
(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
1079.1 1120.7 -529.5 1059.1 464
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 0.08044 0.2836
dadid (Intercept) 0.23517 0.4849
Number of obs: 474, groups: momid, 10; dadid, 15
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.2470 0.2401 -1.029 0.30358
crossStxSp -1.8690 0.3947 -4.736 2.18e-06 ***
crossSpxSp 1.0548 0.3700 2.851 0.00436 **
crossSpxSt 0.2237 0.2588 0.865 0.38729
hr24hr 1.7728 0.1851 9.578 < 2e-16 ***
crossStxSp:hr24hr 0.3882 0.3539 1.097 0.27261
crossSpxSp:hr24hr -0.3402 0.3090 -1.101 0.27085
crossSpxSt:hr24hr -0.4167 0.2610 -1.597 0.11031
---
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 |
|---|---|---|---|---|---|
| SpxSp4hr | SpxSp24hr | -1.433 | 0.247 | 0.000 | *** |
| SpxSt4hr | SpxSt24hr | -1.356 | 0.184 | 0.000 | *** |
| SpxSp4hr | SpxSt4hr | 0.831 | 0.326 | 0.013 | * |
| SpxSp24hr | SpxSt24hr | 0.908 | 0.373 | 0.016 | * |
| StxSp4hr | StxSp24hr | -2.161 | 0.303 | 0.000 | *** |
| StxSt4hr | StxSt24hr | -1.773 | 0.185 | 0.000 | *** |
| StxSp4hr | StxSt4hr | -1.869 | 0.395 | 0.000 | *** |
| StxSp24hr | StxSt24hr | -1.481 | 0.342 | 0.000 | *** |
| SpxSp4hr | StxSt4hr | 1.055 | 0.370 | 0.006 | ** |
| SpxSp24hr | StxSt24hr | 0.715 | 0.410 | 0.081 | . |
| SpxSp4hr + SpxSt4hr | StxSt4hr - StxSp4hr | 3.148 | 0.508 | 0.000 | *** |
| SpxSp24hr + SpxSt24hr | StxSt24hr - StxSp24hr | 2.002 | 0.502 | 0.000 | *** |
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 |
|---|---|---|---|---|
| StxSt | 4hr | 0.44 | 0.50 | 0.38 |
| StxSp | 4hr | 0.11 | 0.15 | 0.08 |
| SpxSp | 4hr | 0.69 | 0.75 | 0.63 |
| SpxSt | 4hr | 0.49 | 0.56 | 0.43 |
| StxSt | 24hr | 0.82 | 0.86 | 0.78 |
| StxSp | 24hr | 0.51 | 0.58 | 0.44 |
| SpxSp | 24hr | 0.90 | 0.93 | 0.87 |
| SpxSt | 24hr | 0.79 | 0.83 | 0.75 |
(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.2-1.
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.2-1},
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},
}