4
$\begingroup$

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?

$\endgroup$
3
  • 1
    $\begingroup$ In R fitting a logistic regression is nothing more than specifying glm(y ~ x1 + x2 + x3 + x4, family=binomial). Does this answer your question? $\endgroup$
    – AdamO
    Feb 18, 2021 at 21:15
  • $\begingroup$ Why are you building your own likelihood function? Why not use R's glm function? $\endgroup$ Feb 18, 2021 at 22:24
  • 1
    $\begingroup$ 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. $\endgroup$
    – mmr09
    Feb 19, 2021 at 0:39

2 Answers 2

5
$\begingroup$

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))
```
$\endgroup$
0
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.