时间戳

首页 > 解决方案 > MLE 错误：“vmmin”中的初始值不是有限的

r - MLE 错误：“vmmin”中的初始值不是有限的

问题描述

我们模拟了一个数据集并创建了一个模型。

set.seed(459)

# seed mass
n <- 1000
seed.mass <- round(rnorm(n, mean = 250, sd = 75),digits = 1)

## Setting up the deterministic function
detFunc <- function(a,b,x){
  return(exp(a+b*x)) / (1+exp(a+b*x))
}

# logit link function for the binomial
inv.link <- function(z){   
  p <-1/(1+exp(-z))
  return(p)
  }

#setting a and b values
a <- -2.109
b <- 0.02

# Simulating data 
germination <- (rbinom(n = n, size = 10,
                        p = inv.link(detFunc(x = seed.mass, a = a, b = b))
                        ))/10

## make data frame
mydata <- data.frame("predictor" = seed.mass, "response" = germination)

# plotting the data
tmp.x <- seq(0,1e3,length.out=500)
plot(germination ~ seed.mass,
     xlab = "seed mass (mg)",
     ylab = "germination proportion")
  lines(tmp.x,inv.link(detFunc(x = tmp.x, a = a, b = b)),col="red",lwd=2)

当我们检查我们创建的模型并推断参数时，我们得到一个错误：

optim(par = c(a = -2.109, b = 0.02), fn = function (p) 中的错误：“vmmin”中的初始值不是有限的

library(bbmle)
  mod1<-mle2(response ~ dbinom(size = 10,
                        p = inv.link(detFunc(x = predictor, a = a, b = b))
                        ),
             data = mydata,
             start = list("a"= -2.109 ,"b"= 0.02))

我们很困惑，无法弄清楚为什么会出现这个错误。

标签： rmodelingmle

解决方案

您的问题是您试图将二项式结果（必须是整数）拟合到一个比例。

您可以将round(response*10)其用作预测器（将比例放回计数范围；round()因为在浮点数学中(a/b)*b并不总是完全等于...）具体来说，使用您的设置a

mod1 <- mle2(round(response*10) ~ dbinom(size = 10,
                    p = inv.link(detFunc(x = predictor, a = a, b = b))
                    ),
         data = mydata,
         start = list(a = -2.109 ,b = 0.02))

工作正常。coef(mod1)是 {-1.85, 0.018}，很可能接近您开始时的真实值（我们不希望准确地恢复真实值，除非作为许多模拟的平均值 [即使这样 MLE 也只是渐近无偏的，即对于大数据集...]

最接近的问题是，尝试dbinom()使用非整数值进行评估会给出NA. 您的模型拟合的完整输出将是：

optim(par = c(a = -2.109, b = 0.02), fn = function (p) 中的错误：'vmmin' 中的初始值不是有限的另外：有 50 个或更多警告（使用 warnings() 来查看前 50 个）

检查那些额外的警告总是一个好主意......在这种情况下，它们都是形式

1：在 dbinom(x = c(1, 1, 1, 0.8, 1, 1, 1, 1, 1, 1, 1, 0.8, ... : 非整数 x = 0.800000

这可能给了你一个线索......

PS你可以使用qlogis()和plogis()从基础R为你的链接和反向链接功能......

推荐阅读