Schiedea salicaria and S. menziesii pollen tube growth
John Powers
1/30/2020
Purpose
Test for post-pollination, pre-zygotic reproductive barriers between two plant species, Schiedea salicaria and S. menziesii by measuring pollen tube growth from parental and hybrid crosses at two time points.
Response variables:
- Counts - assume a zero-inflated quasi-Poisson distribution for overdispersed count data with an excess of zeroes
- 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 beta-binomial distribution for overdispersed 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 beta-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=" ", cache=T)
library(ggplot2)
library(ggbeeswarm)
library(ggsignif)
library(ggforce)
library(tidyr)
library(dplyr)
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, ...)
}Read data
#setwd("~/MyDocs/MEGA/UCI/Schiedea/Analysis/HybridSeeds")
pops.sali <- c("842")
pops.menz <- c("849", "950")
tubes <- read.csv("pollentubes_salmen.csv")
tubes <- within(tubes, {
hr <- relevel(hr,"4 hr")
levels(hr) <- c("4hr","24hr")
ratioraw <- bottom/mid
adjmid <- ifelse(bottom > mid, bottom, mid)
ratio <- bottom/adjmid
justmid <- mid - bottom
pollnum <- factor(pollnum)
mompop <- factor(substr(momid, 1, 3))
dadpop <- factor(substr(dadid, 1, 3))
momcross <- factor(ifelse(mompop %in% pops.sali, "S", "M"))
dadcross <- factor(ifelse(dadpop %in% pops.sali, "S", "M"))
cross <- factor(paste0(momcross,"x",dadcross))
crosshr <- factor(paste0(cross, hr))
})
str(tubes) 'data.frame': 153 obs. of 19 variables:
$ line : int 114 422 423 454 455 330 331 452 453 206 ...
$ pollnum : Factor w/ 153 levels "22","23","24",..: 37 94 95 104 105 76 77 102 103 46 ...
$ cross : Factor w/ 4 levels "MxM","MxS","SxM",..: 4 4 4 4 4 3 3 3 3 4 ...
$ hr : Factor w/ 2 levels "4hr","24hr": 1 1 2 1 2 1 2 1 2 1 ...
$ momid : Factor w/ 18 levels "842 192-37","842 201-46",..: 2 2 2 2 2 2 2 2 2 4 ...
$ dadid : Factor w/ 19 levels "842 192-37","842 201-46",..: 5 5 5 9 9 15 15 18 18 8 ...
$ date : Factor w/ 16 levels "2020-01-07","2020-01-08",..: 2 11 11 12 12 8 8 12 12 4 ...
$ collect_date: Factor w/ 20 levels "","2020-01-07",..: 3 1 1 13 14 9 10 13 14 5 ...
$ mid : int 9 4 9 8 13 6 9 9 12 4 ...
$ bottom : int 3 1 7 3 11 1 9 3 11 1 ...
$ crosshr : Factor w/ 8 levels "MxM24hr","MxM4hr",..: 8 8 7 8 7 6 5 6 5 8 ...
$ dadcross : Factor w/ 2 levels "M","S": 2 2 2 2 2 1 1 1 1 2 ...
$ momcross : Factor w/ 2 levels "M","S": 2 2 2 2 2 2 2 2 2 2 ...
$ dadpop : Factor w/ 3 levels "842","849","950": 1 1 1 1 1 2 2 3 3 1 ...
$ mompop : Factor w/ 3 levels "842","849","950": 1 1 1 1 1 1 1 1 1 1 ...
$ justmid : int 6 3 2 5 2 5 0 6 1 3 ...
$ ratio : num 0.333 0.25 0.778 0.375 0.846 ...
$ adjmid : int 9 4 9 8 13 6 9 9 12 4 ...
$ ratioraw : num 0.333 0.25 0.778 0.375 0.846 ...Inventory of attempted crosses
with(tubes, table(cross, hr)) hr
cross 4hr 24hr
MxM 22 22
MxS 20 20
SxM 18 18
SxS 17 16with(tubes, table(cross, momid)) momid
cross 842 192-37 842 201-46 842 22-4 842 221F 842 37-4 842 43A-41
MxM 0 0 0 0 0 0
MxS 0 0 0 0 0 0
SxM 4 4 4 6 4 4
SxS 4 5 4 4 4 4
momid
cross 842 51-20 842 86-38 842 B99 849 122-3 849 122-3 849 4-17 849 57-39
MxM 0 0 0 4 0 4 4
MxS 0 0 0 4 2 4 4
SxM 4 5 1 0 0 0 0
SxS 4 4 0 0 0 0 0
momid
cross 849 57C 849 61A 849 77-2 849 93C 950 406-2
MxM 6 6 6 8 6
MxS 4 8 4 4 6
SxM 0 0 0 0 0
SxS 0 0 0 0 0with(tubes, table(dadid, momid)) momid
dadid 842 192-37 842 201-46 842 22-4 842 221F 842 37-4 842 43A-41
842 192-37 0 0 0 0 0 0
842 201-46 0 0 0 0 0 2
842 215-2 0 0 0 0 0 0
842 22-4 0 0 0 0 0 0
842 221F 2 3 0 0 0 0
842 251-6 2 0 0 0 0 0
842 37-4 0 0 2 0 0 0
842 43A-41 0 0 0 2 0 0
842 51-20 0 2 2 0 2 2
842 86-38 0 0 0 2 2 0
849 122-3 0 0 2 0 0 0
849 4-17 2 0 0 0 0 2
849 57-39 0 0 0 0 2 0
849 57C 2 0 0 0 0 0
849 61A 0 2 0 2 2 0
849 93C 0 0 2 0 0 2
849122-3 0 0 0 0 0 0
950 406-2 0 2 0 4 0 0
950 501 0 0 0 0 0 0
momid
dadid 842 51-20 842 86-38 842 B99 849 122-3 849 122-3 849 4-17
842 192-37 0 0 0 0 0 0
842 201-46 2 2 0 0 0 0
842 215-2 0 0 0 0 0 0
842 22-4 0 2 0 0 0 0
842 221F 0 0 0 0 0 0
842 251-6 0 0 0 0 0 0
842 37-4 0 0 0 0 0 2
842 43A-41 0 0 0 2 0 0
842 51-20 0 0 0 2 2 2
842 86-38 2 0 0 0 0 0
849 122-3 2 0 0 0 0 0
849 4-17 0 3 0 0 0 0
849 57-39 0 0 0 2 0 2
849 57C 0 0 0 0 0 0
849 61A 2 0 0 2 0 0
849 93C 0 2 1 0 0 0
849122-3 0 0 0 0 0 0
950 406-2 0 0 0 0 0 2
950 501 0 0 0 0 0 0
momid
dadid 849 57-39 849 57C 849 61A 849 77-2 849 93C 950 406-2
842 192-37 0 0 0 2 0 0
842 201-46 2 2 0 0 0 0
842 215-2 0 2 0 0 0 0
842 22-4 0 0 0 0 0 0
842 221F 2 0 4 0 0 2
842 251-6 0 0 0 0 0 0
842 37-4 0 0 2 0 2 0
842 43A-41 0 0 2 2 0 0
842 51-20 0 0 0 0 2 0
842 86-38 0 0 0 0 0 4
849 122-3 0 2 0 0 5 0
849 4-17 0 0 2 0 2 2
849 57-39 0 0 0 2 0 0
849 57C 0 0 0 0 0 0
849 61A 0 2 0 2 0 4
849 93C 2 2 0 2 0 0
849122-3 0 0 0 0 1 0
950 406-2 0 0 4 0 0 0
950 501 2 0 0 0 0 0Plot raw data
Pollen tubes
(tubeplot.mid <-
ggplot(tubes, aes(x=cross, y=mid, color=hr)) +
ggtitle("Pollen tubes at middle of style") +
geom_violin(position=position_dodge(), draw_quantiles=c(0.25,0.5,0.75), size=0.5, bw=0.8) +
geom_point(position=position_jitterdodge(jitter.height=0)) +
scale_color_brewer(palette=3, type="qual") +
scale_y_continuous(limits=c(0,NA)) +
theme_bw() + theme(axis.text.x = element_text(angle=90, hjust = 1, vjust = 0.5)))
(tubeplot.bottom <-
ggplot(tubes, aes(x=cross, y=bottom, color=hr)) +
ggtitle("Pollen tubes at bottom of style") +
geom_violin(position=position_dodge(), draw_quantiles=c(0.25,0.5,0.75), size=0.5, bw=0.8) +
geom_point(position=position_jitterdodge(jitter.height=0)) +
scale_color_brewer(palette=3, type="qual") +
theme_bw() + theme(axis.text.x = element_text(angle=90, hjust = 1, vjust = 0.5)))
#sinaplot
set.seed(3)
(tubeplot.ratio <-
ggplot(tubes, aes(x=cross, y=ratio, color=hr)) +
geom_hline(aes(yintercept=1), size=0.4)+
geom_violin(position=position_dodge(), draw_quantiles=c(0.25,0.5,0.75), size=0.5, bw=0.05) +
geom_sina(bw=0.05)+
labs(x="", y="Ratio of pollen tubes at bottom : middle of style",color="Time") +
scale_color_brewer(palette=3, type="qual", labels=c("4 hr","24 hr"), direction=1) +
scale_y_continuous(expand = expand_scale(add=c(0,0.01)), breaks = scales::pretty_breaks(n = 5)) +
scale_x_discrete(labels=rev(c("S x S", "S x M", "M x S", "M x M"))) +
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), axis.ticks.x = element_blank()))
(tubeplot.justmid <-
ggplot(tubes, aes(x=cross, y=justmid, color=hr)) +
ggtitle("Difference of pollen tubes at middle and bottom of style") +
geom_violin(position=position_dodge(), draw_quantiles=c(0.25,0.5,0.75), size=0.5, bw=0.8) +
geom_point(position=position_jitterdodge(jitter.height=0)) +
ylab("mid - botttom") +
scale_color_brewer(palette=3, type="qual") +
theme_bw() + theme(axis.text.x = element_text(angle=90, hjust = 1, vjust = 0.5)))
(tubeplot.mid.bottom <-
ggplot(tubes, aes(x=bottom, y=mid)) +
ggtitle("Pollen tubes at middle vs. bottom of style") +
geom_abline(slope=1, size=1) +
geom_point(color="grey") +
stat_ellipse(aes(linetype=hr, color=cross), size=1) +
coord_cartesian(xlim= c(0,max(tubes$bottom)), ylim=c(0,max(tubes$mid))) +
theme_bw())
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
709.4 739.7 -344.7 689.4 143
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 5.673e-03 7.532e-02
dadid (Intercept) 1.158e-10 1.076e-05
Number of obs: 153, groups: momid, 18; dadid, 19
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.958601 0.084503 23.178 < 2e-16 ***
crossMxS -0.146921 0.121141 -1.213 0.225
crossSxM -0.029912 0.125946 -0.237 0.812
crossSxS 0.018884 0.126442 0.149 0.881
hr24hr 0.443205 0.102594 4.320 1.56e-05 ***
crossMxS:hr24hr 0.080650 0.153547 0.525 0.599
crossSxM:hr24hr -0.017726 0.154544 -0.115 0.909
crossSxS:hr24hr -0.001807 0.155637 -0.012 0.991
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1signif.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("***", "**", "*", ".", " "))
}
crosshr.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=c("SxS4hr - SxS24hr == 0",
"SxM4hr - SxM24hr == 0",
"SxS4hr - SxM4hr == 0",
"SxS24hr - SxM24hr == 0",
"MxS4hr - MxS24hr == 0",
"MxM4hr - MxM24hr == 0",
"MxS4hr - MxM4hr == 0",
"MxS24hr - MxM24hr == 0",
"SxS4hr - MxM4hr == 0",
"SxS24hr - MxM24hr == 0",
"SxS4hr+SxM4hr - MxM4hr-MxS4hr == 0",
"SxS24hr+SxM24hr - MxM24hr-MxS24hr == 0"))) %>%
summary(test=adjusted(type="fdr")) %>% tidy %>%
separate(lhs, c("lvl1", "lvl2"), sep = " - ", extra="merge") %>%
mutate(signif = p.value < 0.05,
stars = signif.num(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(3, 8,10:12)], digits=c(0,0,3,3,3,8), caption="Contrasts")| lvl1 | lvl2 | estimate | std.error | statistic | p.value | stars |
|---|---|---|---|---|---|---|
| SxS4hr | SxS24hr | -0.441 | 0.117 | -3.771 | 0.00064918 | *** |
| SxM4hr | SxM24hr | -0.425 | 0.116 | -3.681 | 0.00069619 | *** |
| SxS4hr | SxM4hr | 0.049 | 0.127 | 0.383 | 0.88127770 | |
| SxS24hr | SxM24hr | 0.065 | 0.104 | 0.622 | 0.80125766 | |
| MxS4hr | MxS24hr | -0.524 | 0.114 | -4.586 | 0.00005434 | *** |
| MxM4hr | MxM24hr | -0.443 | 0.103 | -4.320 | 0.00009362 | *** |
| MxS4hr | MxM4hr | -0.147 | 0.121 | -1.213 | 0.54048578 | |
| MxS24hr | MxM24hr | -0.066 | 0.095 | -0.697 | 0.80125766 | |
| SxS4hr | MxM4hr | 0.019 | 0.126 | 0.149 | 0.88127770 | |
| SxS24hr | MxM24hr | 0.017 | 0.105 | 0.162 | 0.88127770 | |
| SxS4hr + SxM4hr | MxM4hr - MxS4hr | 0.136 | 0.190 | 0.714 | 0.80125766 | |
| SxS24hr + SxM24hr | MxM24hr - MxS24hr | 0.036 | 0.159 | 0.224 | 0.88127770 |
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 |
|---|---|---|---|---|
| MxM | 4hr | 7.1 | 7.7 | 6.5 |
| MxS | 4hr | 6.1 | 6.7 | 5.6 |
| SxM | 4hr | 6.9 | 7.6 | 6.3 |
| SxS | 4hr | 7.2 | 7.9 | 6.6 |
| MxM | 24hr | 11.0 | 11.8 | 10.3 |
| MxS | 24hr | 10.3 | 11.1 | 9.6 |
| SxM | 24hr | 10.5 | 11.4 | 9.7 |
| SxS | 24hr | 11.2 | 12.2 | 10.4 |
(mid.emplot <- ggplot(cld.mod, aes(y=response, x=cross, fill=hr)) +
geom_signif(y_position=cont$height+13, 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,16), expand = expand_scale(add=c(0,0.5)), breaks = scales::pretty_breaks(n = 8)) +
scale_fill_brewer(palette=3, type="qual", labels=c("24 hr","4 hr"), direction=-1) +
scale_x_discrete(labels=c("S x S", "S x M", "M x S", "M x M")) +
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), axis.ticks.x = element_blank()))
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
618.6 648.9 -299.3 598.6 143
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 5.001e-03 7.072e-02
dadid (Intercept) 1.379e-10 1.174e-05
Number of obs: 153, groups: momid, 18; dadid, 19
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.5190 0.1664 3.119 0.00181 **
crossMxS -0.1113 0.2458 -0.453 0.65063
crossSxM 0.3024 0.2295 1.318 0.18756
crossSxS 0.3363 0.2309 1.457 0.14524
hr24hr 1.7551 0.1780 9.857 < 2e-16 ***
crossMxS:hr24hr 0.0028 0.2660 0.011 0.99160
crossSxM:hr24hr -0.3511 0.2492 -1.409 0.15884
crossSxS:hr24hr -0.3083 0.2509 -1.229 0.21914
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1signif.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("***", "**", "*", ".", " "))
}
crosshr.x <- data.frame(crosshr=levels(tubes$crosshr), x=rep(1:4, each=2)+rep(c(0.25, -0.25), 4))
cont <- glmmTMB(bottom~crosshr + (1|momid) + (1|dadid), data=tubes, family="poisson") %>%
glht(linfct = mcp(crosshr=c("SxS4hr - SxS24hr == 0",
"SxM4hr - SxM24hr == 0",
"SxS4hr - SxM4hr == 0",
"SxS24hr - SxM24hr == 0",
"MxS4hr - MxS24hr == 0",
"MxM4hr - MxM24hr == 0",
"MxS4hr - MxM4hr == 0",
"MxS24hr - MxM24hr == 0",
"SxS4hr - MxM4hr == 0",
"SxS24hr - MxM24hr == 0",
"SxS4hr+SxM4hr - MxM4hr-MxS4hr == 0",
"SxS24hr+SxM24hr - MxM24hr-MxS24hr == 0"))) %>%
summary(test=adjusted(type="fdr")) %>% tidy %>%
separate(lhs, c("lvl1", "lvl2"), sep = " - ", extra="merge") %>%
mutate(signif = p.value < 0.05,
stars = signif.num(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(3, 8,10:12)], digits=c(0,0,3,3,3,8), caption="Contrasts")| lvl1 | lvl2 | estimate | std.error | statistic | p.value | stars |
|---|---|---|---|---|---|---|
| SxS4hr | SxS24hr | -1.447 | 0.177 | -8.182 | 0.00000000 | *** |
| SxM4hr | SxM24hr | -1.404 | 0.174 | -8.051 | 0.00000000 | *** |
| SxS4hr | SxM4hr | 0.034 | 0.222 | 0.153 | 0.87865975 | |
| SxS24hr | SxM24hr | 0.077 | 0.111 | 0.693 | 0.73239841 | |
| MxS4hr | MxS24hr | -1.758 | 0.198 | -8.892 | 0.00000000 | *** |
| MxM4hr | MxM24hr | -1.755 | 0.178 | -9.857 | 0.00000000 | *** |
| MxS4hr | MxM4hr | -0.111 | 0.246 | -0.453 | 0.78075981 | |
| MxS24hr | MxM24hr | -0.109 | 0.102 | -1.059 | 0.49624565 | |
| SxS4hr | MxM4hr | 0.336 | 0.231 | 1.457 | 0.29047577 | |
| SxS24hr | MxM24hr | 0.028 | 0.110 | 0.253 | 0.87264132 | |
| SxS4hr + SxM4hr | MxM4hr - MxS4hr | 0.750 | 0.339 | 2.216 | 0.06411707 | . |
| SxS24hr + SxM24hr | MxM24hr - MxS24hr | 0.088 | 0.166 | 0.528 | 0.78075981 |
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 |
|---|---|---|---|---|
| MxM | 4hr | 1.7 | 2.0 | 1.4 |
| MxS | 4hr | 1.5 | 1.8 | 1.3 |
| SxM | 4hr | 2.3 | 2.7 | 1.9 |
| SxS | 4hr | 2.4 | 2.8 | 2.0 |
| MxM | 24hr | 9.7 | 10.5 | 9.0 |
| MxS | 24hr | 8.7 | 9.4 | 8.1 |
| SxM | 24hr | 9.3 | 10.0 | 8.5 |
| SxS | 24hr | 10.0 | 10.9 | 9.2 |
(bottom.emplot <- ggplot(cld.mod, aes(y=response, x=cross, fill=hr)) +
geom_signif(y_position=cont$height+13, 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,16), expand = expand_scale(add=c(0,0.5)), breaks = scales::pretty_breaks(n = 8)) +
scale_fill_brewer(palette=3, type="qual", labels=c("24 hr","4 hr"), direction=-1) +
scale_x_discrete(labels=c("S x S", "S x M", "M x S", "M x M")) +
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), axis.ticks.x = element_blank()))
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
440.7 471.0 -210.3 420.7 143
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 2.548e-09 5.048e-05
dadid (Intercept) 3.411e-11 5.840e-06
Number of obs: 153, groups: momid, 18; dadid, 19
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.16821 0.18823 -6.206 5.43e-10 ***
crossMxS 0.04761 0.28219 0.169 0.8660
crossSxM 0.46294 0.26808 1.727 0.0842 .
crossSxS 0.43824 0.26922 1.628 0.1036
hr24hr 3.16689 0.27310 11.596 < 2e-16 ***
crossMxS:hr24hr -0.35298 0.39470 -0.894 0.3712
crossSxM:hr24hr -0.47912 0.40061 -1.196 0.2317
crossSxS:hr24hr -0.35748 0.40974 -0.872 0.3830
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1signif.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("***", "**", "*", ".", " "))
}
crosshr.x <- data.frame(crosshr=levels(tubes$crosshr), x=rep(1:4, each=2)+rep(c(0.25, -0.25), 4))
cont <- glmmTMB(cbind(bottom, adjmid-bottom)~crosshr + (1|momid) + (1|dadid), data=tubes, family="binomial") %>%
glht(linfct = mcp(crosshr=c("SxS4hr - SxS24hr == 0",
"SxM4hr - SxM24hr == 0",
"SxS4hr - SxM4hr == 0",
"SxS24hr - SxM24hr == 0",
"MxS4hr - MxS24hr == 0",
"MxM4hr - MxM24hr == 0",
"MxS4hr - MxM4hr == 0",
"MxS24hr - MxM24hr == 0",
"SxS4hr - MxM4hr == 0",
"SxS24hr - MxM24hr == 0",
"SxS4hr+SxM4hr - MxM4hr-MxS4hr == 0",
"SxS24hr+SxM24hr - MxM24hr-MxS24hr == 0"
))) %>%
summary(test=adjusted(type="fdr")) %>% tidy %>%
separate(lhs, c("lvl1", "lvl2"), sep = " - ", extra="merge") %>%
mutate(signif = p.value < 0.05,
stars = signif.num(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(3, 8,10:12)], digits=c(0,0,3,3,3,8), caption="Contrasts")| lvl1 | lvl2 | estimate | std.error | statistic | p.value | stars |
|---|---|---|---|---|---|---|
| SxS4hr | SxS24hr | -2.809 | 0.305 | -9.198 | 0.00000000 | *** |
| SxM4hr | SxM24hr | -2.688 | 0.293 | -9.170 | 0.00000000 | *** |
| SxS4hr | SxM4hr | -0.025 | 0.271 | -0.091 | 0.92742238 | |
| SxS24hr | SxM24hr | 0.097 | 0.325 | 0.298 | 0.92742238 | |
| MxS4hr | MxS24hr | -2.814 | 0.285 | -9.875 | 0.00000000 | *** |
| MxM4hr | MxM24hr | -3.167 | 0.273 | -11.596 | 0.00000000 | *** |
| MxS4hr | MxM4hr | 0.048 | 0.282 | 0.169 | 0.92742238 | |
| MxS24hr | MxM24hr | -0.305 | 0.276 | -1.107 | 0.46026729 | |
| SxS4hr | MxM4hr | 0.438 | 0.269 | 1.628 | 0.20711790 | |
| SxS24hr | MxM24hr | 0.081 | 0.309 | 0.261 | 0.92742238 | |
| SxS4hr + SxM4hr | MxM4hr - MxS4hr | 0.854 | 0.391 | 2.181 | 0.06998126 | . |
| SxS24hr + SxM24hr | MxM24hr - MxS24hr | 0.370 | 0.426 | 0.867 | 0.57852805 |
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 |
|---|---|---|---|---|
| MxM | 4hr | 0.24 | 0.27 | 0.20 |
| MxS | 4hr | 0.25 | 0.29 | 0.21 |
| SxM | 4hr | 0.33 | 0.37 | 0.29 |
| SxS | 4hr | 0.33 | 0.37 | 0.28 |
| MxM | 24hr | 0.88 | 0.90 | 0.86 |
| MxS | 24hr | 0.84 | 0.87 | 0.82 |
| SxM | 24hr | 0.88 | 0.90 | 0.85 |
| SxS | 24hr | 0.89 | 0.91 | 0.86 |
(ratio.emplot <- ggplot(cld.mod, aes(y=response, x=cross, fill=hr)) +
geom_signif(y_position=cont$height/15+0.9, 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 = expand_scale(add=c(0,0.05)), breaks = scales::pretty_breaks(n = 5)) +
scale_fill_brewer(palette=3, type="qual", labels=c("24 hr","4 hr"), direction=-1) +
scale_x_discrete(labels=c("S x S", "S x M", "M x S", "M x M")) +
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), axis.ticks.x = element_blank()))
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 (2019). emmeans: Estimated Marginal Means, aka
Least-Squares Means. R package version 1.4.3.
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 = {2019},
note = {R package version 1.4.3},
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},
}