Note: Change of various wordpress themes and R updates may mean the code below might not work without some minor tweaks. Find the full code here (last checked 20/02/23). 

https://raw.githubusercontent.com/gerard-ricardo/blog/main/Three%20factor%20dose-response%20with%20binomial%20GLMMs%20and%20ECxs.R

This post will describe how to fit binomial GLMMs to a 3 fixed factor x 1 continuous factor design, a common design in wet-lab experiments. This design is my preferred design I used during my Postdoc for multifactor experiments because it accurately predicts the effect size (ECx) while investigating two factors. For example, you may want to know how much stressor X impacts coral survival, and how this changes with a second factor Y.

The easiest way to use this code is to reformat your headings in your data sheet to match the general heading I use here, therefore no code needs to be altered.

Note: This code does not go deep into data exploration, model selection, or model validation to keep things simple. I am going to assume a random ‘tank’ effect (obs), which occurs in almost all experiments.

1. Import some data and label type. Your data should be organised in this type of long form

data1 <- read.table(file=”https://raw.githubusercontent.com/gerard-ricardo/data/master/mockmulticb&#8221;, header= TRUE,dec=”,”, na.strings=c(“”,”.”,”NA”))
head(data1)
options(scipen = 999) # turn off scientific notation

#2)Organising and wrangling
str(data1) #check data type is correct
data1$raw.x <- as.numeric(as.character(data1$raw.x))
# data1$counts <- as.numeric(as.character(data1$counts))
data1$suc <- as.integer(as.character(data1$suc))
data1$tot <- as.integer(as.character(data1$tot))
data1$prop <- data1$suc/data1$tot
data1$obs <- factor(formatC(1:nrow(data1), flag=”0″, width = 3))# unique tank ID for later on
nrow(data1)
str(data1)
data1 = data1[complete.cases(data1), ] #make sure import matches NA type

Okay let’s do some very basic data exploration

table(data1$tot)

> table(data1$tot)

10 12 15 16 17 18 19 20 21 22 23 24 25 26 27 29 30 31 33 34 35 39 40
1 1 3 5 5 5 5 3 5 3 6 4 3 3 4 4 3 3 2 1 1 1 1

This isn’t great, it means each tank has a different number of animals, and predictions for tanks with very few animals will be very coarse. GLMs can weight this unbalance a bit so let’s leave in for the moment. We can see this in the data my making the size of each tank relative to the number of animals in it.

library(ggplot2)
p0 = ggplot()+geom_point(data1, mapping = aes(x = raw.x, y = prop),position = position_jitter(width = .02), alpha = 0.50,size = data1$tot*0.2)+theme_light()
p0 = p0 + facet_wrap(~factor)+scale_x_log10(name =”raw.x”)
p0

Okidok, at this stage you may want to launch into data exploration and consider some of the following

• is the control healthy/ or representative of field data health?
• how does the control relate to the first treatment?
• are data evenly spaced across the logx scale?
• is the effect linear or non-linear?
• is there any non-monotonic (bell shape) effects?
• are there any outliers?

This data set has a number of issues, but lets push on ahead anyway.

#4) Fit the model
library(lme4)
library(splines)
md3 <- glmer(cbind(suc,(tot – suc)) ~ scale(raw.x) * factor + (1|obs) ,family = binomial (link = logit),data = data1)
summary(md3)
library(RVAideMemoire) #GLMM overdispersion test
overdisp.glmer(md3) #Overdispersion for GLMM
#is there an interactive effect?
#is it overdispersed?

At this stage, many would stop and interpret. They would talk about the interactive effect and for significant effects say something like ‘for every unit increase of X, there is a decrease in response of blah blah blah. Neither statement is super useful for Regulators that actually need to make decisions. Why? Because p-values in lab assays often just reflect the experimental design. Experiments that use high concentrations and lots of replicates will likely find significant effects, but these may not be biologically meaningful. And by placing the interpretation on the treatment (which is to be managed), rather than response makes the metric much more usuable. So let’s take an effect-size approach by looking at the magnitude of change of the response. These metrics (ECx) can also be used in other models such as SSDs and meta-analyses. But first let’s predict and plot up the models. Here I’m going to switch to glmmTMB because of convergence issues.

library(glmmTMB)
md1 =glmmTMB(cbind(suc,(tot-suc))~raw.x*factor+(1|obs) ,data1,family=’binomial’ )
summary(md1)
vec.x = seq(min(data1$raw.x), max(data1$raw.x), length = 100)
df.x <- expand.grid(raw.x = vec.x,
factor = levels(data1$factor))
#make sure names are the same as model
mm <- model.matrix(~raw.x*factor, df.x) # build model matrix to sub the parameters into the equation
eta <- mm %*% fixef(md1)$cond #sub in parameters to get mean on logit scale
df.x$prediction <- as.vector(exp(eta) / (1 + exp(eta))) #back-transform mean from logit scale
se <- sqrt(diag(mm %*% vcov(md1)$cond %*% t(mm)))
df.x$upper <- exp(eta + 1.96 *se) /(1 + exp(eta + 1.96 *se)) #work out upper 95% CI
df.x$lower <- exp(eta – 1.96 *se) /(1 + exp(eta – 1.96 *se)) #work out lower 95% CI
data1$factor = factor(data1$factor,levels = c(“a”, “b”, “c”)) #Set levels in order for model
library(ggplot2)
p0= ggplot()
p0= p0+ geom_point(data = data1, aes(x = raw.x, y = prop,alpha = 0.1), color = ‘steelblue’, size = data1$tot*0.1, position=position_jitter(width = .01))
p0= p0+ geom_line(data = df.x, aes(x = raw.x, y = prediction), color = ‘grey30′, size=1)
p0= p0+ geom_ribbon(data = df.x, aes(x = raw.x, ymin=lower, ymax=upper,fill=’grey’), alpha=0.2)
p0= p0+ scale_x_log10(limits = c(0.9, 100))
p0 = p0+ labs(x=expression(Deposited~sediment~(mg~”cm”^{-2})),
y=expression(Recruit~survival~(prop.)))
p0= p0+ scale_y_continuous( limits = c(-0.05, 1.01))
p0= p0+ theme_light()
p0= p0+ scale_fill_manual( values = c(“grey”,”khaki2″))
p0= p0+ theme(legend.position=”none”)
p0= p0+ facet_wrap(~factor, nrow = 1)
p0

You can see that for each factor, there are some thresholds occurring. For the first level of the categorical factor (the left plot) the curve starts to bend steeply around 10 units. For the third plot, it occurs < 10 units, although the confidence bands are still wide in this region.

We want to predict what levels of the continuous stressor (on the x-axis) causes a 50% effect i.e how much of the stressor will impact 50% of the population. We can do this by interpolating onto the curve. I have used Vernables* dose.p function to do this. In this example I want to compare every EC50 against ‘factor a”s EC50 (represented in the red line).

Yep the code below looks intense but let me explain what is fundamentally going on here.

1. We say what type of effect we want to see (here 50% is ec = 50).
2. Work out the top of each model at the control
3. Find 50% from that.
4. Get everything on the logit scale ready for interpolation
5. Interpolate for the mean for each curve. Note the model needs to be releved to get the coefficients for the second curve.
6. Use some variance-co-variance to get the standard errors
7. Use the old Wald approximation to approximate the 95% confidence intervals
8. Wrap everything into a data.frame
9. This needs to be done for each curve, so I ‘relevel’ for each factor.

#6) Getting ECx’s
ec = 50 #put in your ecx here
library(dplyr)
group.fac <-df.x %>% group_by(factor)%>%summarise(estimate = max(prediction))%>%as.data.frame()
top1 = group.fac$estimate[1] #the modelled control of factor 1
inhib1 = top1 -((ec/100) * top1) #an x% decrease from the control for factor 1
library(VGAM)
eta1 <- logit(inhib1)
data1$factor <- relevel(data1$factor, ref = “b”) #set reference levels for GLMs
md2 <- glmmTMB(cbind(suc,(tot-suc))~raw.x*factor+(1|obs) ,data1,family=’binomial’ )
data1$factor <- relevel(data1$factor, ref = “c”) #set reference levels for GLMs
md3 <- glmmTMB(cbind(suc,(tot-suc))~raw.x*factor+(1|obs) ,data1,family=’binomial’ )
data1$factor = factor(data1$factor,levels = c(“a”, “b”, “c”)) #Set levels in order for model
betas1 = fixef(md1)$cond[1:2] #intercept and slope for ref 1
betas2 = fixef(md2)$cond[1:2] #intercept and slope for ref 2
betas3 = fixef(md3)$cond[1:2]
ecx1 <- (eta1 – betas1[1])/betas1[2]
ecx2 <- (eta1 – betas2[1])/betas2[2]
ecx3 <- (eta1 – betas3[1])/betas3[2]
pd1 <- -cbind(1, ecx1)/betas1[2]
pd2 <- -cbind(1, ecx2)/betas2[2]
pd3 <- -cbind(1, ecx3)/betas3[2]
ff1 = as.matrix(vcov(md1)$cond[1:2,1:2])
ff2 = as.matrix(vcov(md2)$cond[1:2,1:2])
ff3 = as.matrix(vcov(md3)$cond[1:2,1:2])
ec.se1 <- sqrt(((pd1 %*% ff1 )* pd1) %*% c(1, 1))
ec.se2 <- sqrt(((pd2 %*% ff2 )* pd2) %*% c(1, 1))
ec.se3 <- sqrt(((pd3 %*% ff3 )* pd3) %*% c(1, 1))
upper1 = (ecx1+ec.se1*1.96)
lower1 = (ecx1-ec.se1*1.96)
upper2 = (ecx2+ec.se2*1.96)
lower2 = (ecx2-ec.se2*1.96)
upper3 = (ecx3+ec.se2*1.96)
lower3 = (ecx3-ec.se2*1.96)
ec.df1 = data.frame(ecx1, lower1, upper1)
ec.df2 = data.frame(ecx2, lower2, upper2)
ec.df3 = data.frame(ecx3, lower3, upper3)
ecall = cbind(ec.df1, ec.df2, ec.df3)
ec.df1 #this is your factor 1 ecx values
ec.df2 #this is your factor 1 ecx values
ec.df3 #
p0= p0+ geom_hline(yintercept = inhib1, col = ‘red’, linetype=”dashed”)
p0

Finally lets visualize the EC50’s on the plot

#Visualising the ECX
ecx.all = bind_cols(data.frame(factor = c(‘a’, ‘b’, ‘c’)), data.frame(ecx = c(ecx1, ecx2, ecx3)), data.frame(inhib = c(inhib1, inhib1, inhib1)))
upper.all = bind_cols(data.frame(factor = c(‘a’, ‘b’, ‘c’)), data.frame(upper = c(upper1, upper2, upper3)), data.frame(inhib = c(inhib1, inhib1, inhib1)))
lower.all = bind_cols(data.frame(factor = c(‘a’, ‘b’, ‘c’)), data.frame(lower = c(lower1, lower2, lower3)), data.frame(inhib = c(inhib1, inhib1, inhib1)))
p0 = p0 + geom_point(data = upper.all, aes(x = upper, y = inhib), color = ‘red’, size=2)
p0 = p0 + geom_point(data = ecx.all, aes(x = ecx, y = inhib), color = ‘red’, size=2)
p0 = p0 + geom_point(data = lower.all, aes(x = lower, y = inhib), color = ‘red’, size=2)
p0

You can see the EC50 values decrease between the factors, especially in factor c.  The EC10 for the left hand plot is quite a bit lower, so a regulator may set lower ‘trigger values’ for areas where these two combinations of stressors occur. The curve on the left has an EC50 of 46.11, but the curve on right has an EC50 > 27.83, so ‘factor c’ is of much greater concern for regulators. In my view this is much better than simply going off p-values. However we can still compare the EC50s statistically, which is a topic for a later post.

I hope this post was useful to you, and will hopefully get researchers new to modelling thinking through the lens of effect size.

*Venables WN, Ripley BD. Modern applied statistics with S-PLUS. Springer Science & Business Media. 2013.

This is a draft bit of code to allows the user to take .csv files that commonly come from image analysis programs, combine them together into a data.frame in R, wrangle the rows and columns into the required format for ordination analysis. This code will also allow for csv files with uneven columns.

Note: This code is in beta testing. Make sure to check your data at each step to make sure it all makes sense.

A Photoquad screenshot

1.  This is a typical csv with taxa in column C

2. Have all your csv’s in a seperate folder

3.  The first chunck of code imports all the csv data into a tibble in R

library(dplyr)
library(readr)
df <- list.files(full.names = TRUE) %>%
lapply(read_csv) %>%
bind_rows
head(df)

4. The second chunk groups the tibble by image and flips to wide format. Cells created by mashing everything together give NA. These can be set to zero as no count was made for that taxa.

df2
library(tidyverse)
df3 = df2 %>% group_by(Image) %>% pivot_wider(names_from = ‘spp Name’, values_from = ‘Cov% per species’, names_prefix = “area”) #year goes to columns, #their areas go as the values, area is the prefix
df3
df3[is.na(df3)] <- 0 #all NA to zero
names(df3)
tail(df3)

5. A basic nMDS from the data

# library(tidyverse)
df4 = subset(df3, select=-Image)
library(vegan)
comMDS2 <- metaMDS(df4, distance = “bray”, autotransform = T) #20 random #interactions
library(ggordiplots)
ggordi <- gg_ordiplot(comMDS2, groups = df3$Image, plot = T) ##

You are probably wondering how to automate the image analysis part. There are quite a few people working on this but also comes with a huge amount of challenges for images high in complexity and diversity. Laboratory (wet -lab ) experimental images are much easier and I will add a post on how to do this at some stage.

After conducting coral spawning experiments since 2013 in wet-labs – a total of 13 coral spawning events – I finally had the opportunity to see it for the first time in the wild. The chance of seeing spawning was uncertain, as water temperatures for the time of year were unseasonably low, and we had heard very few corals had developed pigmented eggs. Still, Magnetic Island corals had spawned without fail in October since I can remember, and so we decided that at worst we would get a few nice days holiday. There wasn’t much to lose.

The first night was really a night designated to test equipment and find a good site, but within a few minutes of decent, we saw the first setting corals. The first coral I observed spawning was at ~8 pm, and various Acropora and Montipora spp. spawned past our ascent time of 9:30pm. There was lots still set as we exited the water.

Note: I haven’t done any editing on the photos so I can keep a reference of how things looked underwater just in case I want to work on some of these spp. in the future. All photos taken with a dive-torch and a TG-6.

The first coral I saw setting around 8pm.

I find this one interesting. The branches have almost no egg-sperm bundles, but are super-dense between the branches.

The coral-eating Drupella snail

The zooplankton attracted to the light were nuts. …and had a good sting.

On the second night there was, suprisingly, no spawning. It still made for a fun night-dive but there was hardly the flurry of activity from the previous night. The plankton, crabs, sharks and fish had mostly disappeared, worn out from the night before.

This blenny apparently likes to sleep above the waterline. 

This text originally appeared in the Annual 2016 Australian Coral Reef Society Newsletter

From the paper: Ricardo, G. F., Negri, A. P., Jones, R. J. & Stocker, R. That sinking feeling: Suspended sediments can prevent the ascent of coral egg bundles. Scientific Reports 6, 21567 (2016). doi:10.1038/srep21567

Fig. 1. Predicted reduction in the number of coral egg-sperm bundles reaching the water surface due to sediment ballasting, and microscopy of the ballasted bundles after failed ascent. A) Percent reduction in bundle ascents. Green and blue contours denote the SS loads that reduce successful ascents by 10% (EC10) and 50% (EC50). B) EC10 values for the ascent failure (EC10, green and purple lines, referring to two different sediment grain radii) and encounter failure (EC10, red line). C) Optical microscopy image showing sediment grains attached to a Montipora digitata bundle. D) Coloured backscatter scanning electron microscopy micrographs of an Acropora nasuta bundle, showing sediment grains in yellow and biotic matter in purple. Scale bar = 200. 

Inshore reefs are regularly exposed to elevated concentrations of suspended sediments associated with natural resuspension events, runoff and dredging. While it has been recognised that suspended sediments can negatively impact coral recruitment, the vulnerability of critical reproductive stages prior to fertilisation have not been considered. In a recent paper, we demonstrate a previously unrecognized pre-fertilization mechanism that can markedly reduce the probability that coral gametes reach the water surface (where fertilisation takes place) in the presence of commonly encountered concentrations of suspended sediments.

Synchronised broadcast spawning of gametes is the dominant mode of reproduction in hard corals and the buoyancy of egg-sperm bundles is critical to maximise gamete encounters and fertilization at the ocean surface. We demonstrate that (i) during their ascent to the surface, the bundles can intercept suspended sediment grains that stick to their mucous coating, and (ii) under a broad range of realistic environmental conditions the ballasting effect of the sediment grains is sufficient to cause a sizeable fraction of bundles to sink and never reach the water surface. The detrimental impact of this loss of ascending bundles on egg-sperm encounters is magnified because the bundles carry both eggs and sperm, and the encounter rate is proportional to the product of their respective concentrations at the surface. Even for bundles that remain positively buoyant, reaching the surface is delayed, further reducing egg-sperm encounter probabilities, and hence fertilisation success. Observations of this mechanism are successfully captured by a mathematical model that predicts the reduction in ascent probability and egg-sperm encounters as a function of sediment load, depth and particle grain size. For reefs at 15 m deep, the model predicts that a coarse silt could reduce 10% of egg-sperm encounters at 35 mg L⁻¹, and for a reef at 5 m deep, reduce 10% of egg-sperm encounters at 106 mg L⁻¹.

This is the first study to examine the effects of environmental pressures on the success of coral gamete ascent, which can have important flow-on effects for the recruitment success of corals following disturbance. This new mechanism intensifies the cumulative risk posed by other sediment and climate stresses on early life history stages of corals and could contribute to recruitment failure on nearby reefs. Furthermore, the mechanism and model we propose provides a blueprint for related processes in other marine organisms, including some echinoderms, molluscs and fish that rely on positively buoyant eggs for fertilization and are thus vulnerable to sediment ballasting.