Schiedea membranacea and S. kauaiensis pollen tube growth
John Powers
2022-05-10
Purpose
Test for post-pollination, pre-zygotic reproductive barriers between two plant species, Schiedea membranacea and S. kauaiensis 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
Excludes rows where there was no pollen on stigmas.
setwd("~/MyDocs/MEGA/UCI/Schiedea/Analysis/HybridSeeds")
mem.sp <- c("mem","864") #dadsp sometimes has 864 instead of mem
tubes <- read_csv("pollentubes_memkau.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") %>%
mutate(status = fct_explicit_na(status, "OK"),
mid = ifelse(status=="no pollen on stigmas", NA, mid), #set to NA to drop, 0 to keep
bottom = ifelse(status=="no pollen on stigmas", NA, bottom)) %>% #see SGW email 2022-05-09, NOTES tab
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,
momcross = ifelse(momsp %in% mem.sp, "M", "K"),
dadcross = ifelse(dadsp %in% mem.sp, "M", "K"),
momid = paste(momcross, momplant),
dadid = paste(dadcross, dadplant),
cross = factor(paste0(momcross,"x",dadcross), levels=c("MxM","MxK","KxK","KxM")),
crosshr = factor(paste0(cross, hr), levels = paste(rep(levels(cross),each=2), levels(hr), sep=""))) %>%
mutate_if(is.character, as.factor)
str(tubes) tibble [621 × 25] (S3: tbl_df/tbl/data.frame)
$ line : num [1:621] 1 3 3 3 3 3 4 4 4 5 ...
$ emasculation_date: Date[1:621], format: "2020-12-05" "2020-12-05" ...
$ pollination_date : Date[1:621], format: "2020-12-08" "2020-12-08" ...
$ momsp : Factor w/ 2 levels "kaua","mem": 2 2 2 2 2 2 2 2 2 2 ...
$ momplant : Factor w/ 12 levels "1ax234-1#7","1ax234-1#9",..: 9 9 9 9 9 9 9 9 9 9 ...
$ dadsp : Factor w/ 3 levels "864","kaua","mem": 1 1 1 1 1 1 1 1 1 2 ...
$ dadplant : Factor w/ 19 levels "12074-7x525-11#2",..: 18 18 18 18 18 18 18 18 18 4 ...
$ vial : num [1:621] 363 365 365 365 365 365 366 366 366 367 ...
$ hr : Factor w/ 2 levels "4hr","24hr": 1 2 2 2 2 2 2 2 2 1 ...
$ count_date : Date[1:621], format: "2021-12-01" "2021-12-01" ...
$ notes : Factor w/ 17 levels "JP: added 0 to bottom_4 per SGW email",..: NA NA NA NA NA NA NA NA NA NA ...
$ status : Factor w/ 3 levels "no pollen on stigmas",..: 3 3 3 3 3 3 3 3 3 3 ...
$ stigma : Factor w/ 5 levels "1","2","3","4",..: 1 1 2 3 4 5 1 2 3 1 ...
$ mid : num [1:621] 1 4 2 3 5 5 5 5 5 0 ...
$ bottom : num [1:621] 0 4 6 5 4 3 4 7 4 0 ...
$ ratioraw : num [1:621] 0 1 3 1.67 0.8 ...
$ adjmid : num [1:621] 1 4 6 5 5 5 5 7 5 0 ...
$ ratio : num [1:621] 0 1 1 1 0.8 0.6 0.8 1 0.8 NaN ...
$ justmid : num [1:621] 1 0 -4 -2 1 2 1 -2 1 0 ...
$ momcross : Factor w/ 2 levels "K","M": 2 2 2 2 2 2 2 2 2 2 ...
$ dadcross : Factor w/ 2 levels "K","M": 2 2 2 2 2 2 2 2 2 1 ...
$ momid : Factor w/ 12 levels "K 1ax234-1#7",..: 10 10 10 10 10 10 10 10 10 10 ...
$ dadid : Factor w/ 19 levels "K 12074-7x525-11#2",..: 19 19 19 19 19 19 19 19 19 4 ...
$ cross : Factor w/ 4 levels "MxM","MxK","KxK",..: 1 1 1 1 1 1 1 1 1 2 ...
$ crosshr : Factor w/ 8 levels "MxM4hr","MxM24hr",..: 1 2 2 2 2 2 2 2 2 3 ...tube.contrasts <- c("KxK4hr - KxK24hr == 0",
"KxM4hr - KxM24hr == 0",
"KxK4hr - KxM4hr == 0",
"KxK24hr - KxM24hr == 0",
"MxK4hr - MxK24hr == 0",
"MxM4hr - MxM24hr == 0",
"MxK4hr - MxM4hr == 0",
"MxK24hr - MxM24hr == 0",
"KxK4hr - MxM4hr == 0",
"KxK24hr - MxM24hr == 0",
"KxK4hr+KxM4hr - MxM4hr-MxK4hr == 0",
"KxK24hr+KxM24hr - MxM24hr-MxK24hr == 0")Inventory of crosses
with(tubes, table(cross, hr)) hr
cross 4hr 24hr
MxM 75 80
MxK 87 95
KxK 73 76
KxM 67 68with(tubes, table(cross, momid)) momid
cross K 1ax234-1#7 K 1ax234-1#9 K 1ax525-11#7 K 2a K 3ax525-1#6 K Prog2x525-1#9
MxM 0 0 0 0 0 0
MxK 0 0 0 0 0 0
KxK 46 12 21 23 24 23
KxM 38 4 24 21 24 24
momid
cross M 2021-1 M 5x9c-7 M 864 192-30 M mem#2x201-11#2 M mem#2x201-11#6
MxM 34 23 37 9 23
MxK 34 31 36 11 37
KxK 0 0 0 0 0
KxM 0 0 0 0 0
momid
cross M mem#2x201-11#7
MxM 29
MxK 33
KxK 0
KxM 0with(tubes, table(dadid, momid)) momid
dadid K 1ax234-1#7 K 1ax234-1#9 K 1ax525-11#7 K 2a K 3ax525-1#6
K 12074-7x525-11#2 0 0 0 0 12
K 1ax234-1#? 0 0 0 0 0
K 1ax234-1#1 0 0 0 0 12
K 1ax234-1#10 0 0 21 0 0
K 1ax234-1#5 21 0 0 0 0
K 1ax234-1#7 0 0 0 11 0
K 1ax234-1#9 0 0 0 0 0
K 1ax525-1#20 0 0 0 0 0
K 1ax525-1#7 0 12 0 0 0
K 2a 0 0 0 0 0
K 3ax525-1#6 0 0 0 0 0
K Prog2x525-1#9 25 0 0 12 0
M 2021-1 19 0 0 9 0
M 5x9c-7 19 0 0 12 0
M 864 192-30 0 0 12 0 12
M mem#2x201-11#1 0 0 0 0 0
M mem#2x201-11#2 0 0 12 0 12
M mem#2x201-11#6 0 0 0 0 0
M mem#2x201-11#7 0 4 0 0 0
momid
dadid K Prog2x525-1#9 M 2021-1 M 5x9c-7 M 864 192-30
K 12074-7x525-11#2 0 0 0 0
K 1ax234-1#? 0 0 14 0
K 1ax234-1#1 0 0 0 0
K 1ax234-1#10 0 0 0 0
K 1ax234-1#5 0 0 0 0
K 1ax234-1#7 11 0 0 0
K 1ax234-1#9 0 0 0 0
K 1ax525-1#20 12 0 0 0
K 1ax525-1#7 0 0 0 0
K 2a 0 17 0 0
K 3ax525-1#6 0 0 0 18
K Prog2x525-1#9 0 17 17 18
M 2021-1 12 0 15 0
M 5x9c-7 12 17 0 0
M 864 192-30 0 17 0 0
M mem#2x201-11#1 0 0 0 0
M mem#2x201-11#2 0 0 0 37
M mem#2x201-11#6 0 0 8 0
M mem#2x201-11#7 0 0 0 0
momid
dadid M mem#2x201-11#2 M mem#2x201-11#6 M mem#2x201-11#7
K 12074-7x525-11#2 0 0 0
K 1ax234-1#? 0 0 0
K 1ax234-1#1 0 0 0
K 1ax234-1#10 11 0 17
K 1ax234-1#5 0 0 0
K 1ax234-1#7 0 20 0
K 1ax234-1#9 0 0 16
K 1ax525-1#20 0 0 0
K 1ax525-1#7 0 0 0
K 2a 0 0 0
K 3ax525-1#6 0 0 0
K Prog2x525-1#9 0 17 0
M 2021-1 0 19 0
M 5x9c-7 0 4 0
M 864 192-30 0 0 0
M mem#2x201-11#1 0 0 10
M mem#2x201-11#2 0 0 19
M mem#2x201-11#6 0 0 0
M mem#2x201-11#7 9 0 0tubes %>% pivot_longer(c(momid, dadid), values_to = "plant", values_ptypes = character()) %>%
count(plant, name) %>% pivot_wider(values_from="n") %>% arrange(plant) %>%
kable(caption="Number of stigmas split by paternal or maternal plant")| plant | dadid | momid |
|---|---|---|
| K 12074-7x525-11#2 | 12 | |
| K 1ax234-1#? | 14 | |
| K 1ax234-1#1 | 12 | |
| K 1ax234-1#10 | 49 | |
| K 1ax234-1#5 | 21 | |
| K 1ax234-1#7 | 42 | 84 |
| K 1ax234-1#9 | 16 | 16 |
| K 1ax525-1#20 | 12 | |
| K 1ax525-1#7 | 12 | |
| K 1ax525-11#7 | 45 | |
| K 2a | 17 | 44 |
| K 3ax525-1#6 | 18 | 48 |
| K Prog2x525-1#9 | 106 | 47 |
| M 2021-1 | 74 | 68 |
| M 5x9c-7 | 64 | 54 |
| M 864 192-30 | 41 | 73 |
| M mem#2x201-11#1 | 10 | |
| M mem#2x201-11#2 | 80 | 20 |
| M mem#2x201-11#6 | 8 | 60 |
| M mem#2x201-11#7 | 13 | 62 |
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
2451.3 2495.6 -1215.7 2431.3 611
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 0.1290 0.3591
dadid (Intercept) 0.0837 0.2893
Number of obs: 621, groups: momid, 12; dadid, 19
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.64237 0.20278 3.168 0.00154 **
crossMxK 0.14063 0.18357 0.766 0.44363
crossKxK 0.06958 0.27734 0.251 0.80192
crossKxM -1.66026 0.29197 -5.686 1.30e-08 ***
hr24hr 0.89526 0.09306 9.620 < 2e-16 ***
crossMxK:hr24hr -0.55386 0.13162 -4.208 2.58e-05 ***
crossKxK:hr24hr -0.03299 0.13108 -0.252 0.80132
crossKxM:hr24hr 1.10261 0.21535 5.120 3.06e-07 ***
---
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 |
|---|---|---|---|---|---|
| KxK4hr | KxK24hr | -0.862 | 0.092 | 0.000 | *** |
| KxM4hr | KxM24hr | -1.998 | 0.194 | 0.000 | *** |
| KxK4hr | KxM4hr | 1.730 | 0.250 | 0.000 | *** |
| KxK24hr | KxM24hr | 0.594 | 0.173 | 0.001 | ** |
| MxK4hr | MxK24hr | -0.341 | 0.093 | 0.001 | *** |
| MxM4hr | MxM24hr | -0.895 | 0.093 | 0.000 | *** |
| MxK4hr | MxM4hr | 0.141 | 0.184 | 0.591 | |
| MxK24hr | MxM24hr | -0.413 | 0.168 | 0.021 | * |
| KxK4hr | MxM4hr | 0.070 | 0.277 | 0.880 | |
| KxK24hr | MxM24hr | 0.037 | 0.264 | 0.890 | |
| KxK4hr + KxM4hr | MxM4hr - MxK4hr | -1.731 | 0.482 | 0.001 | *** |
| KxK24hr + KxM24hr | MxM24hr - MxK24hr | -0.108 | 0.441 | 0.880 |
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 |
|---|---|---|---|---|
| MxM | 4hr | 1.9 | 2.3 | 1.6 |
| MxK | 4hr | 2.2 | 2.6 | 1.8 |
| KxK | 4hr | 2.0 | 2.5 | 1.7 |
| KxM | 4hr | 0.4 | 0.5 | 0.3 |
| MxM | 24hr | 4.7 | 5.6 | 3.8 |
| MxK | 24hr | 3.1 | 3.7 | 2.6 |
| KxK | 24hr | 4.8 | 5.8 | 4.0 |
| KxM | 24hr | 2.7 | 3.3 | 2.2 |
(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
1769.1 1813.4 -874.5 1749.1 611
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 0.1104 0.3323
dadid (Intercept) 0.1914 0.4375
Number of obs: 621, groups: momid, 12; dadid, 19
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.81464 0.27377 -2.976 0.002924 **
crossMxK 0.24478 0.30857 0.793 0.427614
crossKxK 0.17920 0.36915 0.485 0.627364
crossKxM -1.95535 0.51619 -3.788 0.000152 ***
hr24hr 2.08424 0.16963 12.287 < 2e-16 ***
crossMxK:hr24hr -1.44586 0.23965 -6.033 1.61e-09 ***
crossKxK:hr24hr -0.02679 0.23719 -0.113 0.910078
crossKxM:hr24hr 0.63468 0.49153 1.291 0.196627
---
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 |
|---|---|---|---|---|---|
| KxK4hr | KxK24hr | -2.057 | 0.166 | 0.000 | *** |
| KxM4hr | KxM24hr | -2.719 | 0.461 | 0.000 | *** |
| KxK4hr | KxM4hr | 2.135 | 0.526 | 0.000 | *** |
| KxK24hr | KxM24hr | 1.473 | 0.261 | 0.000 | *** |
| MxK4hr | MxK24hr | -0.638 | 0.169 | 0.000 | *** |
| MxM4hr | MxM24hr | -2.084 | 0.170 | 0.000 | *** |
| MxK4hr | MxM4hr | 0.245 | 0.309 | 0.570 | |
| MxK24hr | MxM24hr | -1.201 | 0.254 | 0.000 | *** |
| KxK4hr | MxM4hr | 0.179 | 0.369 | 0.684 | |
| KxK24hr | MxM24hr | 0.152 | 0.307 | 0.684 | |
| KxK4hr + KxM4hr | MxM4hr - MxK4hr | -2.021 | 0.659 | 0.003 | ** |
| KxK24hr + KxM24hr | MxM24hr - MxK24hr | 0.033 | 0.439 | 0.940 |
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 | 0.4 | 0.6 | 0.3 |
| MxK | 4hr | 0.6 | 0.7 | 0.4 |
| KxK | 4hr | 0.5 | 0.7 | 0.4 |
| KxM | 4hr | 0.1 | 0.1 | 0.0 |
| MxM | 24hr | 3.6 | 4.5 | 2.8 |
| MxK | 24hr | 1.1 | 1.3 | 0.9 |
| KxK | 24hr | 4.1 | 5.1 | 3.4 |
| KxM | 24hr | 1.0 | 1.2 | 0.7 |
(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
1430.9 1475.2 -705.4 1410.9 611
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
momid (Intercept) 0.05553 0.2356
dadid (Intercept) 0.26840 0.5181
Number of obs: 621, groups: momid, 12; dadid, 19
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.16538 0.30272 -3.850 0.000118 ***
crossMxK 0.13916 0.37143 0.375 0.707914
crossKxK -0.02293 0.40595 -0.056 0.954955
crossKxM -0.57841 0.58504 -0.989 0.322823
hr24hr 2.28199 0.22281 10.242 < 2e-16 ***
crossMxK:hr24hr -1.83835 0.31042 -5.922 3.18e-09 ***
crossKxK:hr24hr 0.20474 0.31437 0.651 0.514862
crossKxM:hr24hr -1.04762 0.61086 -1.715 0.086350 .
---
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 |
|---|---|---|---|---|---|
| KxK4hr | KxK24hr | -2.487 | 0.222 | 0.000 | *** |
| KxM4hr | KxM24hr | -1.234 | 0.570 | 0.073 | . |
| KxK4hr | KxM4hr | 0.556 | 0.626 | 0.563 | |
| KxK24hr | KxM24hr | 1.808 | 0.338 | 0.000 | *** |
| MxK4hr | MxK24hr | -0.444 | 0.220 | 0.087 | . |
| MxM4hr | MxM24hr | -2.282 | 0.223 | 0.000 | *** |
| MxK4hr | MxM4hr | 0.139 | 0.371 | 0.772 | |
| MxK24hr | MxM24hr | -1.699 | 0.323 | 0.000 | *** |
| KxK4hr | MxM4hr | -0.023 | 0.406 | 0.955 | |
| KxK24hr | MxM24hr | 0.182 | 0.346 | 0.720 | |
| KxK4hr + KxM4hr | MxM4hr - MxK4hr | -0.740 | 0.689 | 0.484 | |
| KxK24hr + KxM24hr | MxM24hr - MxK24hr | 0.255 | 0.416 | 0.720 |
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.30 | 0.19 |
| MxK | 4hr | 0.26 | 0.32 | 0.22 |
| KxK | 4hr | 0.23 | 0.29 | 0.19 |
| KxM | 4hr | 0.15 | 0.24 | 0.09 |
| MxM | 24hr | 0.75 | 0.80 | 0.70 |
| MxK | 24hr | 0.36 | 0.41 | 0.31 |
| KxK | 24hr | 0.79 | 0.82 | 0.75 |
| KxM | 24hr | 0.38 | 0.45 | 0.31 |
(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},
}