dbinom - 4 parameter logistic regression

Question

library(grid)
library(gridExtra)
library(broom)
library(BiodiversityR)
library("vegan")#[1]
library("MASS")#[2]
library(nlme)#[3]
library("bbmle")

Here are the data

I'm assessing which model best suit my data (null model/glm-poisson/4 parameter log). The idea of using log is to detect at which point the response (number, proportion of species) decreases/increases at certain values of forest cover in the landscape. I've been using the next code to fit a four parameter logistic regression using dpois (y=count of species):

logip=function(p,lambda,x){
  a=p[1]
  b=p[2]
  c=p[3]
  d=p[4]
  Riq1 = d+(a/(1+exp((b-(FOREST700+km))/c)))
  -sum(dpois(x,lambda=Riq1, log=TRUE))
}
parnames(logip)=c("a","b","c","d")

modTR.log=mle2(minuslog=logip, start= c(a=2,b=60,c=3,d=0.1),
               data=list(x=Patch_Richness))

But now I want to use the same approach for dependent variable that is a proportion (y=proportion of species registered in one site). I think I should be using binomial, so I'm trying dbinom in my previous function:

logip=function(p,size,prob){
 a=p[1]
 b=p[2]
 c=p[3]
 d=p[4]
 Riq1 = d+(a/(1+exp((b-(FOREST500+km))/c)))
 -sum(dbinom(size,prob=Riq1))
}parnames(logip)=c("a","b","c","d")

modTR.log=mle2(minuslog=logip, start= c(a=1,b=72,c=3,d=0.1),
               data=list(x=cbind(Regional_Richness,Patch_Richness)))

I'm getting this message:

Error in mle2(minuslog = logip, start = c(a = 1, b = 72, c = 3, d = 0.1), : some named arguments in 'start' are not arguments to the specified log-likelihood function.

I don´t know if it is correct to use dbinom and how to apply this in the function I'm using. Hope you can help me.

For starters you wouldn't use a logistic regression to model Poisson count data. Logistic regression models binary dependent variables that take values 0 or 1, whereas Poisson distributed variables can take on any non-negative integer value. It would be helpful if you could include some of your data, what you're trying to accomplish, and what packages you're using. — LMc
is there a typo in the second function; you use dnorm instead of dbinom — user20650
Google sheets is asking me to ask for access. Can you make the data freely downloadable? — Ben Bolker

Jellyfish Jellyfish · Accepted Answer · 2021-03-05T04:41:55

For the binomial distribution, there are two parameters. It looks like the Patch_Richness is never greater than 2, so I set the size parameter to 2, and used your formula to predict the probability parameter. Notice the log likelihood is -23.

library(bbmle)
text="Bioma_MAPBIOMAS   km  Regional_Richness   Patch_Richness  Richness_prop   FOREST500
Cerrado 35.1    2   2   1   100
Cerrado 131.4   2   2   1   100
Cerrado 40  2   1   0.5 100
Cerrado 8   1   1   1   72.37
Cerrado 28  1   0   0   85.06
Cerrado 5   1   0   0   29.65
Cerrado 5   1   0   0   25.38
Cerrado 28  1   0   0   77.97
Cerrado 5   1   0   0   70.09
Cerrado 28  1   0   0   100
Cerrado 20  1   0   0   97.48
Cerrado 8   1   0   0   66.89
Cerrado 8   1   0   0   77.96
Cerrado 8   1   0   0   65.17
Cerrado 8   1   0   0   50.86
Cerrado 20  1   0   0   89.1
Cerrado 3   1   1   1   31.49
Cerrado 27.8    1   1   1   62.9"
df=read.table(text=text, header=TRUE, stringsAsFactors = FALSE)
logip=function(p,x){
  a=p[[1]]
  b=p[[2]]
  c=p[[3]]
  d=p[[4]]
  Riq1 = d+a/(1+exp((b-(x$FOREST500+x$km))/c))
  if (any(Riq1>=1) | any(Riq1<=0)) {
    return(9999999)
  }
  -sum(log(dbinom(x$Patch_Richness, 2, prob=exp(Riq1)/(1+exp(Riq1)))))
}
parnames(logip)=c("a","b","c","d")

modTR.log=mle2(minuslog=logip, start= c(a=.1,b=72,c=1,d=.1),
               data=list(x=df))
Call:
mle2(minuslogl = logip, start = c(a = 0.1, b = 72, c = 1, d = 0.1), 
    data = list(x = df))

Coefficients:
           a            b            c            d 
3.811005e-02 7.200033e+01 1.000584e+00 8.823313e-04 

Log-likelihood: -22.45

Here is the same thing for the poisson distribution. Note the log likelihood is -14. Thus, the poisson distribution is better given your equation Riq1 = d+(a/(1+exp((b-(x$FOREST500+x$km))/c))) and initial conditions.

logip=function(p,x){
  a=p[[1]]
  b=p[[2]]
  c=p[[3]]
  d=p[[4]]
  Riq1 = d+(a/(1+exp((b-(x$FOREST500+x$km))/c)))
  if (any(Riq1<=0)) {
    return(9999999)
  }
  -sum(dpois(x$Patch_Richness,lambda=Riq1, log=TRUE))
}
parnames(logip)=c("a","b","c","d")

modTR.log=mle2(minuslog=logip, start= c(a=2,b=60,c=3,d=0.1),
               data=list(x=df))
modTR.log

Call:
mle2(minuslogl = logip, start = c(a = 2, b = 60, c = 3, d = 0.1), 
    data = list(x = df))

Coefficients:
          a           b           c           d 
 0.49350452 77.68600468  0.04856004  0.14285921 

Log-likelihood: -14.5

dbinom - 4 parameter logistic regression

1 Answers