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?

• In R fitting a logistic regression is nothing more than specifying glm(y ~ x1 + x2 + x3 + x4, family=binomial). Does this answer your question? Feb 18, 2021 at 21:15
• Why are you building your own likelihood function? Why not use R's glm function? Feb 18, 2021 at 22:24
• Not really. I'm looking for this: F(x) = D+(A-D)/(1+(x/C)^B) where: A = Minimum asymptote. B = Hill's slope. The Hill's slope refers to the steepness of the curve. It could either be positive or negative. C = Inflection point. The inflection point is defined as the point on the curve where the curvature changes direction or signs. D = Maximum asymptote. Feb 19, 2021 at 0:39

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\left(\frac{y}{1-y}\right) = 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))
$$$$


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
`