how to fit four parameter logistic regression (poisson) R? I'm trying to fit a four parameter logistic regression to model bird species richness (Patch_Richness) in response to forest cover (FOREST500). I need to add km as a co-variable to the model (km= kilometers), I'm trying the following (code) but not sure if "km" is properly included in the formula, I just added km after the predictor (FOREST500)
logip=function(p,lambda,x){
  a=p[1]
  b=p[2]
  c=p[3]
  d=p[4]
  Riq1 = d+(a/(1+exp((b-(FOREST500+km))/c)))
  -sum(dpois(x,lambda=Riq1, log=TRUE))
}
parnames(logip)=c("a","b","c","d")
modTR.log=mle2(minuslog=logip, start= c(a=30,b=30, c=3,d=20), data=list(x=Patch_Richness))
summary(modTR.log)

Any suggestions or comments? Is this ok?
 A: I assume that you'd like to fit a four parameter logistic model extended to multiple independent variables.  I think this essentially means that you'd like to do a logistic regression with a floor and ceiling.
4PL:
$$y = d + \frac{a-d}{1 + (\frac{x}{c})^b}$$
In this context, each of the parameters has a specific interpretation.  It doesn't exactly work when you add more independent variables.  I couldn't find any specific references on how to extend the 4PL, so this is how I would accomplish your goal:
Typical Logistic Regression:
$$ln(\frac{y}{1-y}) = X\beta$$
$$y = \frac{e^{X\beta}}{1+e^{X\beta}}$$
Logistic Regression with a floor and ceiling:
$$y = a + (d - a)\frac{e^{X\beta}}{1+e^{X\beta}}$$
You can fit this with a non-linear least squares approach:
set.seed(25401)

# Simulate Data
f <- function(Forest500, km, a, d)
{
  a + (d - a) * exp(-1+2*Forest500+3*km) / (1 + exp(-1 + 2*Forest500 + 3*km))
}

X <- expand.grid(Forest500 = seq(-10, 10, length = 100), km = seq(-10, 10, length = 25))
Patch_Richness <- f(X$Forest500, X$km, 0.2, 0.75) + rnorm(nrow(X), 0, 0.01)

min(Patch_Richness)
max(Patch_Richness)

plot(X$Forest500, Patch_Richness)
plot(X$km, Patch_Richness)

nls(Patch_Richness ~ a + (d - a) * exp(beta0 + beta1 * Forest500 + beta2 * km) / (1 + exp(beta0 + beta1 * Forest500 + beta2 * km)), data = X,
    start = list(a = 0.1, d = 0.9, beta0 = 0.1, beta1 = 0.1, beta2 = 0.1))
```

A: Here's an example comparing drc::drm and gnlm::gnlr:
library(drc)
spinach1 <- subset(spinach, CURVE==1)
model.drm <- drm(SLOPE~DOSE, CURVE,
              fct=LL.4(names=c("B", "D", "A","C")),data=spinach1)
summary(model.drm)


library(gnlm)
attach(spinach1)
model.gnlr <- gnlr(y = SLOPE, 
                   mu =~ (A-D)/(1+(DOSE/C)^B) + D,
                   pmu = list(A=0.1,B=0.1,C=0.1,D=0.1),
                   pshape=log(0.05))
model.gnlr

