时间戳

首页 > 解决方案 > dbinom - 4 参数逻辑回归

r - dbinom - 4 参数逻辑回归

问题描述

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

这是数据

我正在评估哪种模型最适合我的数据（空模型/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以及如何在我正在使用的功能中应用它。希望您能够帮助我。

标签： rfunctionlogistic-regression

解决方案

对于二项分布，有两个参数。看起来 Patch_Richness 永远不会大于 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

推荐阅读