What is pmodels in drc package

Question

Sorry if it's a dumb question, but I am having trouble figuring out how to use pmodels in the drc package. I've searched everywhere online and all I can find is the definition, which is: "a data frame with a many columns as there are parameters in the non-linear function. Or a list containing a formula for each parameter in the nonlinear function." There are examples online, but I have no what it represents. For example, for the commands:

sel.m2 <- drm(dead/total~conc, type, weights=total, data=selenium, fct=LL.2(),
    type="binomial", pmodels=list(~1, ~factor(type)-1))

met.as.m1<-drm(gain ~ dose, product, data = methionine, fct = AR.3(),
    pmodels = list(~1, ~factor(product), ~factor(product)))
    plot(met.as.m1, log = "", ylim = c(1450, 1800))

auxins.m1 <- boxcox(drm(y ~ dose, h, pmodels = data.frame(h, h, 1, h), fct = LL.4(), data = auxins), method = "anova")

I see pmodels as a list and data frame, but what does the "-1"vs "~1" mean or what does it mean to list a factor, what's the significance of the order within the parenthesis?

Answer 1

I agree that it's not well explained for new people. Unfortunately, I can only answer you in part. A late response but for anyone else:

Two resources are available for reference with drc: a) The writers published about drc. See main text and supplementary (S3 in this example) DOI:10.1371/journal.pone.0146021 b) See the drc.pdf and ctrl+f for pmodel to inspect the various uses.

data.frame vs. list depends on the grouping level I believe.

After playing around with my data (subsets), I found that pmodels() = parameter/pooled models aka how you set those parameters to equal (i.e., global/shared or not). With your last example using the auxins df

library(drc)
auxins.m1 <- boxcox(drm(y ~ dose, h, pmodels = data.frame(h, h, 1, h), 
fct = LL.4(), data = auxins), method = "anova")

## changed names to familiar terms by a non-statistician 
auxins.m1 <- boxcox(drm(y ~ dose, h, pmodels = data.frame(h, h, 1, h), 
fct = LL.4(names=c("hill.slope","bot","top","ed50"), data = auxins), method = "anova")

Shows that the top is set to 1. The order is the same as the LL.4(names...)

So if you set

pmodels = data.frame(h, 1, 1, h)     ## ("hill.slope","bot","top","ed50")

as they do in the drc.pdf on pg.10, you'll see that it's to set a common/shared bottom and top.

Check out pg.9 of their supplementary article, it shows that for LL.2, the two-parameter logistic fit has pre-set top = 1 and bottom = 0. The output of

selenium.LL.2.2 <- drm(dead/total~conc, type, weights = total,
data = selenium, fct = LL.2(), type="binomial",
pmodels = list(~factor(type)-1, ~1))       ## ("hill-slope", "ed50")

Shows that ed50 is assumed constant. Alternatively from pg.91 of the drc.pdf:

## Fitting the model with freely varying ED50 values
mecter.free <- drm(rgr ~ dose, pct, data = mecter,
fct = LL.4(), pmodels = list(~1, ~1, ~1, ~factor(pct) - 1))

Unfortunately, it's really not clear what the object-1 means vs. just the object. A better approach might be to use the base drm() without the special case of LL.#()

Check

getMeanFunctions()

to see all available functions

if you're trying to fix a value at a certain value you can

fct = LL.4(fixed = c(NA,0,1,NA)) 
## effectively becomes the standard LL.2()

## or
fct = LL.4(fixed = c(1,0,NA,NA)) 
## common hill slope = 1; assumes baseline correction hence = 0

Related in part; see a lot of drm functions laid out: https://stackoverflow.com/a/39257095