r语言 - dbinom - 4参数逻辑回归

Question

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

这里是数据

我正在评估哪个模型最适合我的数据(null模型/glm-poisson/4参数日志)。使用原木的想法是检测在景观中特定森林覆盖值的响应(物种数量、比例)在哪个点减少/增加。我一直在使用下一个代码来拟合使用dpois (y=物种计数)的四参数逻辑回归:

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))

但是现在我想对因变量使用相同的方法，这是一个比例(y=一个站点中注册的物种的比例)。我想我应该使用二项式，所以我在之前的函数中尝试dbinom:

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)))

我收到这个消息:

mle2(minuslog = logip, start = c(a = 1, b = 72, c = 3, d = 0.1)中的错误:'start'中的一些命名参数不是指定的对数似然函数的参数。

我不知道使用dbinom是否正确，以及如何将其应用于我正在使用的功能。希望你能帮助我。

Answer 1

对于二项分布，有两个参数。看起来patch_abundance永远不会大于2，所以我将size参数设置为2，并使用你的公式来预测概率参数。注意，日志概率是-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 

泊松分布也是如此。注意，对数似然是-14。因此，在方程Riq1 = d+(a/(1+exp((b-(x$FOREST500+x$km))/c)))和初始条件下，泊松分布较好。

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