Consider adding an argument to fix parameter values in models

[darrenkoppel]

Many concentration-response curves I generate use data normalised to a control response. This gives response values in percent with controls defined as being 100%. As such it would be useful to fix a model's upper asymptote to 100% (or 1). This would be equivalent to the 'fixed' argument in the drc package.

Thanks!

[beckyfisher]

Hi Darren, can you provide a reproducible example of doing this in drc?

[darrenkoppel]

Hi Becky,

Sure thing, here's an example using the earthworms dataset from drc package:

library(drc)
data<-earthworms
str(data)

These data are from an avoidance bioassay - two containers are set next to each other, one containing contaminated soil and the other uncontaminated soil. If the earthworms don't like the contaminated soil they will move to the uncontaminated soil.
The set-up is replicated with increasing doses of a contaminant to give a concentration - response gradient. The number of earthworms in the contaminated soil at the end of the bioassay (number) is compared to the number at the start (total) to give an earthworm proportion in the contaminated soil. At dose 0 the expected response should be 0.5 (no avoidance) between contaminated and uncontaminated soils.

We can fit a a 3-parameter log-logistic model:

model_1 <- drm(number/total~dose, weights = total, data = data,
                     fct = LL.3(), type = "binomial")
plot(model_1)

The upper limit/asymptote is a parameter estimated by the model, but we know a priori that at dose 0 the number/total should be 0.5. So we can use the 'fixed' argument of drm to set that model parameter to 0.5:

model_2 <- drm(number/total~dose, weights = total, data = data,
               fct = LL.3(fixed=c(NA,0.5,NA)), type = "binomial")
plot(model_2)

[beckyfisher]

Thanks Darren

[beckyfisher]

Hi Darren, apologies Diego and I have been too busy to implement this as yet. However, I have worked out some code to show you how to use bayesnec to fit the desired model, pull the brms formula and priors, and then refit using brms, using the constant. You can do this with:

library(drc)
data<-earthworms
str(data)

fit <- bnec(number | trials(total) ~ crf(dose, "ecxll3"),
               data, iter = 2e3, chains=1)
priors_fit <- prior_summary(fit$fit, all = FALSE)
priors_fit

inspect priors and use as a guide to build custom priors with constant 0.5 for top

my_prior <- c(prior_string("constant(0.5)", class = "b", nlpar = "top"),
              prior_string("gamma(5, 2.6)", nlpar = "ec50", lb=0.019, ub = 6.11),
              prior_string("normal(0, 5)", nlpar = "beta"))

get the brmsfit

fit_brms <- pull_brmsfit(fit)

update with custom priors

fit2 <- update(fit_brms, prior = my_prior)

note top is now fixed at 0.5

summary(fit2)