Schiedea kaalae and S. hookeri pollen tube growth
John Powers
2022-05-11
Purpose
Test for post-pollination, pre-zygotic reproductive barriers between two plant species, Schiedea kaalae and S. hookeri 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 genotype, specified by its cross and plant number
- dadid - paternal genotype, specified by its cross and plant number
Explore data
library(knitr)
knitr::opts_chunk$set(comment=" ", warning=F, cache=T)
options(knitr.kable.NA = '')
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")
tubes <- read_csv("pollentubes_kaahoo_extra.csv") %>% mutate(across(ends_with("sp"), 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, bottom) %>% #one entry has mid but no bottom count
mutate(vial=as.character(vial),
momcross = toupper(str_sub(momsp,1,1)),
dadcross = toupper(str_sub(dadsp,1,1)),
observer="SGW") %>%
bind_rows(read_csv("pollentubes_kaahoo.csv") %>% #edited from pollentubes.csv
separate(cross, into=c("momcross","dadcross"), sep="x") %>%
mutate(observer="SZ")) %>%
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(momcross, momid),
dadid = paste(dadcross, dadid),
cross = factor(paste0(momcross,"x",dadcross), levels=c("KxK","KxH","HxH","HxK")),
crosshr = factor(paste0(cross, hr), levels = paste(rep(levels(cross),each=2), levels(hr), sep=""))) %>%
mutate_if(is.character, as.factor)
str(tubes) tibble [515 × 24] (S3: tbl_df/tbl/data.frame)
$ line : num [1:515] 1 1 1 2 2 2 3 3 4 4 ...
$ emasculation_date: Date[1:515], format: "2022-01-07" "2022-01-07" ...
$ pollination_date : Date[1:515], format: "2022-01-10" "2022-01-10" ...
$ momsp : Factor w/ 1 level "hook": 1 1 1 1 1 1 1 1 1 1 ...
$ momid : Factor w/ 12 levels "H 794-2","H 866-1",..: 2 2 2 2 2 2 2 2 2 2 ...
$ dadsp : Factor w/ 2 levels "hook","kaal": 1 1 1 1 1 1 1 1 1 1 ...
$ dadid : Factor w/ 20 levels "H 3--15","H 794-2",..: 2 2 2 2 2 2 2 2 2 2 ...
$ vial : Factor w/ 415 levels "100A","100B",..: 38 38 38 42 42 42 46 46 50 50 ...
$ hr : Factor w/ 2 levels "4hr","24hr": 1 1 1 1 1 1 2 2 2 2 ...
$ count_date : Date[1:515], format: "2022-03-01" "2022-03-01" ...
$ notes : Factor w/ 3 levels "JP: changed cross from KxH > KxK",..: NA NA NA NA NA NA NA NA NA NA ...
$ stigma : Factor w/ 5 levels "1","2","3","4",..: 1 2 3 1 2 3 1 2 1 2 ...
$ mid : num [1:515] 0 0 2 0 0 0 4 4 5 5 ...
$ bottom : num [1:515] 0 0 1 0 0 0 3 5 9 4 ...
$ momcross : Factor w/ 2 levels "H","K": 1 1 1 1 1 1 1 1 1 1 ...
$ dadcross : Factor w/ 2 levels "H","K": 1 1 1 1 1 1 1 1 1 1 ...
$ observer : Factor w/ 2 levels "SGW","SZ": 1 1 1 1 1 1 1 1 1 1 ...
$ grains : num [1:515] NA NA NA NA NA NA NA NA NA NA ...
$ ratioraw : num [1:515] NaN NaN 0.5 NaN NaN NaN 0.75 1.25 1.8 0.8 ...
$ adjmid : num [1:515] 0 0 2 0 0 0 4 5 9 5 ...
$ ratio : num [1:515] NaN NaN 0.5 NaN NaN NaN 0.75 1 1 0.8 ...
$ justmid : num [1:515] 0 0 1 0 0 0 1 -1 -4 1 ...
$ cross : Factor w/ 4 levels "KxK","KxH","HxH",..: 3 3 3 3 3 3 3 3 3 3 ...
$ crosshr : Factor w/ 8 levels "KxK4hr","KxK24hr",..: 5 5 5 5 5 5 6 6 6 6 ...tube.contrasts <- c("HxH4hr - HxH24hr == 0",
"HxK4hr - HxK24hr == 0",
"HxH4hr - HxK4hr == 0",
"HxH24hr - HxK24hr == 0",
"KxH4hr - KxH24hr == 0",
"KxK4hr - KxK24hr == 0",
"KxH4hr - KxK4hr == 0",
"KxH24hr - KxK24hr == 0",
"HxH4hr - KxK4hr == 0",
"HxH24hr - KxK24hr == 0",
"HxH4hr+HxK4hr - KxK4hr-KxH4hr == 0",
"HxH24hr+HxK24hr - KxK24hr-KxH24hr == 0")Inventory of crosses
with(tubes, table(cross, hr, observer)) , , observer = SGW
hr
cross 4hr 24hr
KxK 0 0
KxH 0 0
HxH 36 35
HxK 38 39
, , observer = SZ
hr
cross 4hr 24hr
KxK 62 55
KxH 63 38
HxH 46 26
HxK 43 34with(tubes, table(cross, momid, observer)) , , observer = SGW
momid
cross H 794-2 H 866-1 H 866-2 H 879 2-2 H 879 N-5 H 899 2E-1 K 116-2 K 116-3
KxK 0 0 0 0 0 0 0 0
KxH 0 0 0 0 0 0 0 0
HxH 25 22 0 0 0 24 0 0
HxK 27 24 0 0 0 26 0 0
momid
cross K 140-6 K 892-3 K 892-4 K 892-5
KxK 0 0 0 0
KxH 0 0 0 0
HxH 0 0 0 0
HxK 0 0 0 0
, , observer = SZ
momid
cross H 794-2 H 866-1 H 866-2 H 879 2-2 H 879 N-5 H 899 2E-1 K 116-2 K 116-3
KxK 0 0 0 0 0 0 13 0
KxH 0 0 0 0 0 0 33 21
HxH 0 0 11 61 0 0 0 0
HxK 0 0 17 39 21 0 0 0
momid
cross K 140-6 K 892-3 K 892-4 K 892-5
KxK 0 24 56 24
KxH 21 23 0 3
HxH 0 0 0 0
HxK 0 0 0 0with(tubes, table(dadid, momid)) momid
dadid H 794-2 H 866-1 H 866-2 H 879 2-2 H 879 N-5 H 899 2E-1 K 116-2
H 3--15 0 0 0 16 0 0 0
H 794-2 0 10 0 0 0 12 0
H 866-2 12 0 0 45 0 12 0
H 879 2-2 0 0 0 0 0 0 9
H 879 G 2-2 0 0 11 0 0 0 0
H 879 N-5 0 0 0 0 0 0 14
H 892-1 0 0 0 0 0 0 0
H 892-3 0 0 0 0 0 0 10
H 899 2E-1 13 12 0 0 0 0 0
K 115-2 0 0 0 0 0 0 0
K 116-1 13 12 0 0 0 0 0
K 116-2 0 12 0 0 0 0 0
K 116-3 0 0 0 0 0 14 0
K 140-6 0 0 3 5 0 0 0
K 169-7 0 0 0 10 0 0 0
K 3587-1 0 0 0 0 0 0 13
K 879-2 0 0 0 15 18 0 0
K 892-3 0 0 0 6 3 0 0
K 892-4 0 0 0 0 0 12 0
K 892-5 14 0 14 3 0 0 0
momid
dadid K 116-3 K 140-6 K 892-3 K 892-4 K 892-5
H 3--15 0 0 0 0 0
H 794-2 0 0 0 0 0
H 866-2 0 0 23 0 3
H 879 2-2 0 0 0 0 0
H 879 G 2-2 0 0 0 0 0
H 879 N-5 20 21 0 0 0
H 892-1 1 0 0 0 0
H 892-3 0 0 0 0 0
H 899 2E-1 0 0 0 0 0
K 115-2 0 0 0 36 0
K 116-1 0 0 0 0 0
K 116-2 0 0 24 0 0
K 116-3 0 0 0 14 24
K 140-6 0 0 0 6 0
K 169-7 0 0 0 0 0
K 3587-1 0 0 0 0 0
K 879-2 0 0 0 0 0
K 892-3 0 0 0 0 0
K 892-4 0 0 0 0 0
K 892-5 0 0 0 0 0tubes %>% pivot_longer(c(momid, dadid), values_to = "genotype", values_ptypes = character()) %>%
count(genotype, name) %>% pivot_wider(values_from="n") %>% arrange(genotype) %>%
kable(caption="Number of stigmas split by paternal or maternal genotype")| genotype | dadid | momid |
|---|---|---|
| H 3–15 | 16 | |
| H 794-2 | 22 | 52 |
| H 866-1 | 46 | |
| H 866-2 | 95 | 28 |
| H 879 2-2 | 9 | 100 |
| H 879 G 2-2 | 11 | |
| H 879 N-5 | 55 | 21 |
| H 892-1 | 1 | |
| H 892-3 | 10 | |
| H 899 2E-1 | 25 | 50 |
| K 115-2 | 36 | |
| K 116-1 | 25 | |
| K 116-2 | 36 | 46 |
| K 116-3 | 52 | 21 |
| K 140-6 | 14 | 21 |
| K 169-7 | 10 | |
| K 3587-1 | 13 | |
| K 879-2 | 33 | |
| K 892-3 | 9 | 47 |
| K 892-4 | 12 | 56 |
| K 892-5 | 31 | 27 |
Plot 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
2356.1 2398.5 -1168.0 2336.1 505
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 0.3110 0.5577
dadid (Intercept) 0.1884 0.4341
Number of obs: 515, groups: momid, 12; dadid, 20
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.45293 0.29999 4.843 1.28e-06 ***
crossKxH -0.06835 0.26864 -0.254 0.79915
crossHxH -0.14991 0.44239 -0.339 0.73471
crossHxK -1.18796 0.38278 -3.103 0.00191 **
hr24hr 0.12162 0.07449 1.633 0.10254
crossKxH:hr24hr -0.03452 0.11629 -0.297 0.76658
crossHxH:hr24hr 0.40353 0.09862 4.092 4.28e-05 ***
crossHxK:hr24hr 1.15566 0.14760 7.830 4.89e-15 ***
---
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),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 |
|---|---|---|---|---|---|
| HxH4hr | HxH24hr | -0.525 | 0.065 | 0.000 | *** |
| HxK4hr | HxK24hr | -1.277 | 0.127 | 0.000 | *** |
| HxH4hr | HxK4hr | 1.038 | 0.248 | 0.000 | *** |
| HxH24hr | HxK24hr | 0.286 | 0.240 | 0.466 | |
| KxH4hr | KxH24hr | -0.087 | 0.089 | 0.565 | |
| KxK4hr | KxK24hr | -0.122 | 0.074 | 0.246 | |
| KxH4hr | KxK4hr | -0.068 | 0.269 | 0.799 | |
| KxH24hr | KxK24hr | -0.103 | 0.264 | 0.799 | |
| HxH4hr | KxK4hr | -0.150 | 0.442 | 0.799 | |
| HxH24hr | KxK24hr | 0.254 | 0.441 | 0.799 | |
| HxH4hr + HxK4hr | KxK4hr - KxH4hr | -1.270 | 0.715 | 0.227 | |
| HxH24hr + HxK24hr | KxK24hr - KxH24hr | 0.324 | 0.706 | 0.799 |
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 |
|---|---|---|---|---|
| KxK | 4hr | 4.3 | 5.8 | 3.2 |
| KxH | 4hr | 4.0 | 5.3 | 3.0 |
| HxH | 4hr | 3.7 | 5.0 | 2.7 |
| HxK | 4hr | 1.3 | 1.7 | 1.0 |
| KxK | 24hr | 4.8 | 6.5 | 3.6 |
| KxH | 24hr | 4.4 | 5.8 | 3.2 |
| HxH | 24hr | 6.2 | 8.4 | 4.6 |
| HxK | 24hr | 4.7 | 6.1 | 3.6 |
(mid.emplot <- ggplot(cld.mod, aes(y=response, x=cross, fill=hr)) +
geom_signif(y_position=cont$height+max(cld.mod$uSE), 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,NA), 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
1986.6 2029.1 -983.3 1966.6 505
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 0.4547 0.6743
dadid (Intercept) 0.5061 0.7114
Number of obs: 515, groups: momid, 12; dadid, 20
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.2431 0.4056 3.065 0.002175 **
crossKxH -1.5379 0.4365 -3.523 0.000426 ***
crossHxH -0.4829 0.6069 -0.796 0.426201
crossHxK -2.9025 0.5257 -5.521 3.36e-08 ***
hr24hr 0.3079 0.0842 3.657 0.000255 ***
crossKxH:hr24hr 0.6037 0.1784 3.385 0.000713 ***
crossHxH:hr24hr 0.9775 0.1234 7.922 2.34e-15 ***
crossHxK:hr24hr 1.9668 0.2669 7.368 1.73e-13 ***
---
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),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 |
|---|---|---|---|---|---|
| HxH4hr | HxH24hr | -1.285 | 0.090 | 0.000 | *** |
| HxK4hr | HxK24hr | -2.275 | 0.253 | 0.000 | *** |
| HxH4hr | HxK4hr | 2.420 | 0.427 | 0.000 | *** |
| HxH24hr | HxK24hr | 1.430 | 0.385 | 0.000 | *** |
| KxH4hr | KxH24hr | -0.912 | 0.157 | 0.000 | *** |
| KxK4hr | KxK24hr | -0.308 | 0.084 | 0.001 | *** |
| KxH4hr | KxK4hr | -1.538 | 0.436 | 0.001 | *** |
| KxH24hr | KxK24hr | -0.934 | 0.429 | 0.044 | * |
| HxH4hr | KxK4hr | -0.483 | 0.607 | 0.465 | |
| HxH24hr | KxK24hr | 0.495 | 0.604 | 0.465 | |
| HxH4hr + HxK4hr | KxK4hr - KxH4hr | -1.848 | 0.933 | 0.064 | . |
| HxH24hr + HxK24hr | KxK24hr - KxH24hr | 0.493 | 0.893 | 0.581 |
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 |
|---|---|---|---|---|
| KxK | 4hr | 3.5 | 5.2 | 2.3 |
| KxH | 4hr | 0.7 | 1.1 | 0.5 |
| HxH | 4hr | 2.1 | 3.2 | 1.4 |
| HxK | 4hr | 0.2 | 0.3 | 0.1 |
| KxK | 24hr | 4.7 | 7.1 | 3.1 |
| KxH | 24hr | 1.9 | 2.8 | 1.2 |
| HxH | 24hr | 7.7 | 11.7 | 5.1 |
| HxK | 24hr | 1.9 | 2.7 | 1.3 |
(bottom.emplot <- ggplot(cld.mod, aes(y=response, x=cross, fill=hr)) +
geom_signif(y_position=cont$height+max(cld.mod$uSE), 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,NA), 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
1637.5 1680.0 -808.8 1617.5 505
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 0.1652 0.4065
dadid (Intercept) 0.4298 0.6556
Number of obs: 515, groups: momid, 12; dadid, 20
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.9038 0.3395 2.662 0.00776 **
crossKxH -2.2759 0.4466 -5.096 3.47e-07 ***
crossHxH -0.8376 0.4972 -1.685 0.09205 .
crossHxK -2.7801 0.4470 -6.219 5.00e-10 ***
hr24hr 0.8585 0.1858 4.621 3.82e-06 ***
crossKxH:hr24hr 0.4220 0.2766 1.526 0.12710
crossHxH:hr24hr 1.0905 0.2386 4.571 4.86e-06 ***
crossHxK:hr24hr 0.6268 0.3416 1.835 0.06654 .
---
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 |
|---|---|---|---|---|---|
| HxH4hr | HxH24hr | -1.949 | 0.151 | 0.000 | *** |
| HxK4hr | HxK24hr | -1.485 | 0.287 | 0.000 | *** |
| HxH4hr | HxK4hr | 1.942 | 0.435 | 0.000 | *** |
| HxH24hr | HxK24hr | 2.406 | 0.400 | 0.000 | *** |
| KxH4hr | KxH24hr | -1.280 | 0.206 | 0.000 | *** |
| KxK4hr | KxK24hr | -0.859 | 0.186 | 0.000 | *** |
| KxH4hr | KxK4hr | -2.276 | 0.447 | 0.000 | *** |
| KxH24hr | KxK24hr | -1.854 | 0.444 | 0.000 | *** |
| HxH4hr | KxK4hr | -0.838 | 0.497 | 0.110 | |
| HxH24hr | KxK24hr | 0.253 | 0.506 | 0.673 | |
| HxH4hr + HxK4hr | KxK4hr - KxH4hr | -1.342 | 0.709 | 0.078 | . |
| HxH24hr + HxK24hr | KxK24hr - KxH24hr | -0.046 | 0.652 | 0.943 |
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 |
|---|---|---|---|---|
| KxK | 4hr | 0.71 | 0.78 | 0.64 |
| KxH | 4hr | 0.20 | 0.27 | 0.15 |
| HxH | 4hr | 0.52 | 0.60 | 0.43 |
| HxK | 4hr | 0.13 | 0.18 | 0.10 |
| KxK | 24hr | 0.85 | 0.89 | 0.80 |
| KxH | 24hr | 0.48 | 0.57 | 0.39 |
| HxH | 24hr | 0.88 | 0.91 | 0.84 |
| HxK | 24hr | 0.40 | 0.48 | 0.33 |
(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:
Lenth R (2022). _emmeans: Estimated Marginal Means, aka Least-Squares
Means_. R package version 1.7.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 V. Lenth},
year = {2022},
note = {R package version 1.7.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},
}