Two factor dose-response with binomial GLMMs and ECx’s

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).

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=”;, header= TRUE,dec=”,”, na.strings=c(“”,”.”,”NA”))
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
data1 = data1[complete.cases(data1), ] #make sure import matches NA type

Okay let’s do some very basic data exploration


> 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.

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”)

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
md3 <- glmer(cbind(suc,(tot – suc)) ~ scale(raw.x) * factor + (1|obs) ,family = binomial (link = logit),data = data1)
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.

md1 =glmmTMB(cbind(suc,(tot-suc))~raw.x*factor+(1|obs) ,data1,family=’binomial’ )
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
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})),
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)


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
group.fac <-df.x %>% group_by(factor)%>%summarise(estimate = max(prediction))%>
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
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”)

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)


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.


Automated image point-count to ordination

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

df <- list.files(full.names = TRUE) %>%
lapply(read_csv) %>%

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.

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[] <- 0 #all NA to zero

5. A basic nMDS from the data

# library(tidyverse)
df4 = subset(df3, select=-Image)
comMDS2 <- metaMDS(df4, distance = “bray”, autotransform = T) #20 random #interactions
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.

Coral Spawning 2019

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. 


Suspended sediments can ballast coral egg-sperm bundles

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

Sperm limitation Paper All figures copy

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−1.

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.


Sediments stick and sink coral sperm

Ricardo, G. F., Jones, R. J., Clode, P. L., Humanes, A. & Negri, A. P. Suspended sediments limit coral sperm availability. Scientific Reports (2015).

I began my PhD in 2013 with my research project examining how dredge sediments might impact coral spawning and early developmental stages of corals. Generally, the research project was divided into three parts; how sediment impacts i) fertilisation, ii) larval development and iii) settlement. Coral spawning only occurs once or twice a year on a given reef, but represents an important opportunity for the population to replenish or recover from disturbance events like cyclones, and so management has usually taken a precautious approach by preventing turbidity-generating practises during these rare events.

Sperm limitation Paper All figures copy

Many different pathways in which sediment could impact coral fertilisation

The first question we wanted to help resolve was why previous research had sometimes found sediments were associated with low success rates of fertilisation, but other times not. The variation in fertilisation responses makes it difficult to determine what might be a safe (or risky) sediment concentration while the corals are spawning. We figured that, if we could understand how sediments actually led to a decrease in fertilisation, then we would have some insight into why there had been a range of responses, and perhaps even what sediments presented the greatest risk to coral fertilisation.

1-eggs and embryos (2)

Coral eggs and embryos

So first we designed an experiment based on an ecotoxicology assay by (Marshall 2006) to help us determine whether sediment was impacting the sperm, or the eggs of the coral. We found that while the eggs were capable of fertilising in the presence of high sediment concentrations, far more sperm were needed to achieve an adequate level of fertilisation success. We then conducted sperm counts at the water surface (where the eggs float), which revealed a decrease in sperm numbers in presence of sediment. This indicated that the sediment was actually preventing the sperm from reaching the egg.

ImageJ=1.49oimages=50 slices=50 fps=25 loop=false

Video-microscopy still of coral sperm congregating around a sediment floc

The other thing we discovered when we first ran this experiment was that small flocs appeared on the bottom of the fertilisation containers just after we added the sperm. This gave us a pretty good hunch that the sediments might be sticking to the sperm and cause them to sink away from the eggs. We then used a range of microscopy techniques including scanning electron microscopy, video-microscopy and more recently fluorescence microscopy  to understand interactions between the sediment and the coral gametes. Through a process of confirmation (and elimination) we were able to confidently determine that sperm flocing and sinking was the dominant mechanism affecting coral fertilisation.

Our results have a few implications. Our results suggest that sediments may shrink the coral fertilisation window – which is a brief 1–2-hour period when sperm and eggs can fertilise before wind and waves spread them apart. The coral fertilisation window is nonlinear, meaning that while gametes remain in a saturated concentration, fertilisation remains high. However, when the gamete concentration passes under a certain threshold, fertilisation rates decrease sharply. Sediment can further decrease the sperm concentration. At high sperm concentrations this probably won’t be a problem, but if the sperm concentration is relatively low, the sediments are likely to have a strong impact on fertilisation rates.

Also, because the fertilisation response observed during sediment exposure is highly dependent on the sperm concentration, we urge that a more standardised approach to what sperm concentration is used in lab assays is needed. A lower sperm concentration will ultimately lead to a stronger response, whereas a higher concentration will buffer against the impact of the sediment. Preferably a single standard concentration or a range of concentrations should be used in absence of not knowing the in situ (in field) sperm concentrations.

Orpheus sperm21 (2)

We are finding some sediment types cause floccing more than others

We hope that this research will open up a few avenues for future research questions. It is very likely other marine organisms that broadcast spawn (some fish, echinoderms, molluscs) may also suffer from sediment-sperm floccing. One interesting finding was that carbonate sediment typical of most offshore reefs had no discernible impact on fertilisation rates and we are now investigating which types of sediments cause sperm to floc out of suspension, and therefore present the greatest risk to coral spawning events.

This work was funded by the Western Australian Marine Science Institution and partners.


Marshall DJ (2006) Reliably estimating the effect of toxicants on fertilization success in marine broadcast spawners. Mar Pollut Bull 52:734-738


Masters by Research

Masters by Research – University of Wollongong

Meroplankton larval release and supply in temperate saltmarsh and mangrove habitats

I decided while overseas that I wanted to try my hand at research. At this stage I’d hoped to keep my career options open between terrestrial science and marine science – so I chose a happy medium in estuarine ecology. I was also interested in applied science, and wanted involve some of the amazing marine parks I had access to, with some of my field sites back home in Jervis Bay.

Washing the plankton in the flooded mangrove forest

My main supervisors Todd Minchinton and Andy Davis advised me little was known about zooplankton patterns in Australian estuarine habitats. While trying to develop a research question, I learnt that many invertebrates time their reproduction to environmental cues. I’d heard about this early with the well-known example of multispecific coral spawning in the tropics (more on this later) and also with some fish. I immediately took an interest in this research question as the stakes were quite high for my newfound muddy friends. Because the saltmarsh is only flooded for a brief period every two weeks, their larval release would need to be timed perfectly to allow their young to export to the sea and avoid being stranded in evaporating pools. The surge of larvae at certain times could be also be an important source of food to small fish – making an energetic link between the plants and the larger fish. A simple energetic pathway could be i) saltmarsh photosynthesises solar energy to sugars ii) saltmarsh decomposes into detritus, iii) adult crabs forage on detritus and release larvae iii) small fish eat crab larvae, iv) medium-sized fish eat small fish and iv) large fish eats medium sized fish v) humans eat large fish with chips at Huskisson wharf.

Preparing mentally for night-sampling

The flooded mangrove forests of temperate regions are beautiful. In the tropics, they are littered with prop roots of the Red Mangrove in addition to a constant threat of crocs; but in the south the Grey Mangrove dominates, providing a thick canopy and easy passage for a young plankton tower. The main environmental cues I was looking at were sunlight, moonlight, tidal direction and tidal size. Sampling during the tidal cycles was time-dependant but relatively straight forward; however, sampling at night was quite…interesting and in some cases a little spooky. Luckily I had a great intern that kept me company for most of the sampling period.

I never got sick of looking at the weird creatures one sees under the microscope. Most invisible to the eye, animal plankton show intricate structures at 20 times magnification. But finally after 1.5 years of sampling and counting the larvae, I finally finished the data collection part of the project. I decided two years later to write it up into a manuscript and it was published in the journal Marine Biology