0
votes

We simulated a data set and created a model.

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)

When we check the model we created and infer the parameters, we get an error:

Error in optim(par = c(a = -2.109, b = 0.02), fn = function (p) : initial value in 'vmmin' is not finite

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

We're stumped and can't figure out why we're getting this error.

1
I haven't chased this down fully and won't be able to tonight, but my guess about what's happening is that mle2 is adding log=TRUE to your dbinom() which is making some of the values -Inf and if you look through the code, it's actually calculating minuslogl which would produce a value of Inf. I added log=TRUE into the function call above and got an error that log matched multiple arguments which corroborates the idea that they're adding that into the dbinom() function.DaveArmstrong

1 Answers

1
votes

Your problem is that you're trying to fit a binomial outcome (which must be an integer) to a proportion.

You can use round(response*10) as your predictor (to put the proportion back on the count scale; round() is because (a/b)*b is not always exactly equal to a in floating-point math ...) Specifically, with your setup

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

works fine. coef(mod1) is {-1.85, 0.018}, plausibly close to the true values you started with (we don't expect to recover the true values exactly, except as the average of many simulations [and even then MLE is only asymptotically unbiased, i.e. for large data sets ...]

The proximal problem is that trying to evaluate dbinom() with a non-integer value gives NA. The full output from your model fit would have been:

Error in optim(par = c(a = -2.109, b = 0.02), fn = function (p) : initial value in 'vmmin' is not finite In addition: There were 50 or more warnings (use warnings() to see the first 50)

It's always a good idea to check those additional warnings ... in this case they are all of the form

1: In dbinom(x = c(1, 1, 1, 0.8, 1, 1, 1, 1, 1, 1, 1, 0.8, ... : non-integer x = 0.800000

which might have given you a clue ...

PS you can use qlogis() and plogis() from base R for your link and inverse-link functions ...