mle2 start parameters don't match lower and upper bounds (using L-BFGS-B)

[MPG-GP]

When you pass lower and upper bounds to mle2 for use with L-BFGS-B, the bounds do not automatically line up with the formula and start values.

Here is some code to illustrate the problem.

# Create test data frame with 14 species, 2 canopy levels, and Height measurements every 0.1 meter
spp_levels <- c("abba","acru","acsa","acsp","alin","amel","beal","bepa","coco","cose","frni","pist","poba","potr")
canopy_levels <- c("Open","Shaded")
test_data <- expand.grid(Spp = spp_levels , CanopyClass = canopy_levels , Height = seq(0.1, 20, length.out = 20/0.1) , stringsAsFactors = FALSE)
test_data$Spp <- as.factor(test_data$Spp)
test_data$CanopyClass <- as.factor(test_data$CanopyClass )

# Set test parameters
Tmax_t <- 20
Hmax_t <- 2.5
b_t <- 0.6
overdisp_t <- 0.6

# Create complex power Ricker function with a negative binomial distribution:
# power ricker function:
test_data$twigs <- Tmax_t * Hmax_t^(-b_t*Hmax_t) * exp(b_t*Hmax_t) * test_data$Height^(b_t*Hmax_t) * exp(-b_t*test_data$Height) 
# Add negative binomial random noise:
test_data$twigs_r <- rnbinom( n = length(test_data$twigs), mu = test_data$twigs, size = overdisp_t )

plot(twigs_r ~ Height, data = test_data)

# Use mle2 to fit a model:
mod_test <- mle2( twigs_r ~ dnbinom(mu = Tmax * Hmax^(-b*Hmax) * exp(b*Hmax) * Height^(b*Hmax) * exp(-b*Height), 
           size = overdisp ), data = test_data,  
	     parameters = list(Tmax ~ 1, Hmax ~ 1 , b ~ 1 , overdisp ~ 1 ), 
           optimizer = "optimx", method="L-BFGS-B", start = list( Tmax = 17 , Hmax = 0.9, b = 0.5, overdisp = 2.2),
	     lower = c(Tmax=0.1, Hmax = 0.5, b = 0.01, overdisp = 0.01),
		 upper =  c(Tmax=200, Hmax = 10, b = 3.0 , overdisp = 10.0),
		 control=list( trace = 3, maxit = 100, factr = 1e12), itnmax =100  )

# mod_test converges just fine.

Now, if I try to use mle2 to fit more complex model:

mod_test2 <- mle2( twigs_r ~ dnbinom(mu = Tmax * Hmax^(-b*Hmax) * exp(b*Hmax) * Height^(b*Hmax) * exp(-b*Height), 
           size = overdisp ), data = test_data,  
	     parameters = list(Tmax ~ Spp:CanopyClass -1, Hmax ~ Spp -1 , b ~ Spp -1 , overdisp ~ Spp -1 ), 
           optimizer = "optimx", method="L-BFGS-B", start = list( Tmax = 17 , Hmax = 0.9, b = 0.5, overdisp = 2.2),
	     lower = c(Tmax=0.1, Hmax = 0.5, b = 0.01, overdisp = 0.01),
		 upper =  c(Tmax=200, Hmax = 10, b = 3.0 , overdisp = 10.0),
		 control=list( trace = 3, maxit = 100, factr = 1e12), itnmax =100  )

If throws an error:

fn is  fn1 
gr NULL, npar > 50, kkt set FALSE
Looking for method =  L-BFGS-B 
Methods to be used:[1] "L-BFGS-B"
Function has  70  arguments
par[ 1 ]:  0.1   <? 17   <? 200     In Bounds  
par[ 2 ]:  0.5   <? 17   <? 10     Out of Bounds HIGH  
par[ 3 ]:  0.01   <? 17   <? 3     Out of Bounds HIGH  
par[ 4 ]:  0.01   <? 17   <? 10     Out of Bounds HIGH  
Error in if ((lower[i] <= par[i]) && (par[i] <= upper[i])) { : 

However, if I manually set up the bounds to match the formula parameters:

mod_test3 <- mle2( twigs_r ~ dnbinom(mu = Tmax * Hmax^(-b*Hmax) * exp(b*Hmax) * Height^(b*Hmax) * exp(-b*Height), 
           size = overdisp ), data = test_data,  
	     parameters = list(Tmax ~ Spp:CanopyClass -1, Hmax ~ Spp -1 , b ~ Spp -1 , overdisp ~ Spp -1 ), 
           optimizer = "optimx", method="L-BFGS-B", start = list( Tmax = 17 , Hmax = 0.9, b = 0.5, overdisp = 2.2),
	     lower = c(Tmax=rep(0.1,28) , Hmax = rep(0.5,14) , b = rep(0.01,14), overdisp = rep(0.01,14)),
		 upper =  c(Tmax=rep(200,28), Hmax = rep(10,14), b = rep(3.0,14) , overdisp = rep(10.0,14)),
		 control=list( trace = 3, maxit = 400, factr = 1e12), itnmax =100  )

It runs just fine (although it takes a little while).

But it gets weirder.

If I try to run a model with Hmax ~ 1 (no species variation), I get an even weirder mismatch between start parameters and bounds, even if the bounds are "theoretically" aligned with the formula and the start list.

mod_test4 <- mle2( twigs_r ~ dnbinom(mu = Tmax * Hmax^(-b*Hmax) * exp(b*Hmax) * Height^(b*Hmax) * exp(-b*Height), 
           size = overdisp ), data = test_data,  
	     parameters = list(Tmax ~ Spp:CanopyClass -1, Hmax ~ 1 , b ~ Spp -1 , overdisp ~ Spp -1 ), 
           optimizer = "optimx", method="L-BFGS-B", start = list( Tmax = 17 , Hmax = 0.9, b = 0.5, overdisp = 2.2),
	     lower = c(Tmax=rep(0.1,28) , Hmax = rep(0.5,1) , b = rep(0.01,14), overdisp = rep(0.01,14)),
		 upper =  c(Tmax=rep(200,28), Hmax = rep(10,1), b = rep(3.0,14) , overdisp = rep(10.0,14)),
		 control=list( trace = 3, maxit = 400, factr = 1e12), itnmax =100  )

This code gives an extensive error.

Finally, I found that if you try to manually align the parameters, that doesn't always work, either. For example:

mod_test5 <- mle2( twigs_r ~ dnbinom(mu = Tmax * Hmax^(-b*Hmax) * exp(b*Hmax) * Height^(b*Hmax) * exp(-b*Height), 
           size = overdisp ), data = test_data,  
	     parameters = list(Tmax ~ Spp:CanopyClass -1, Hmax ~ 1 , b ~ Spp -1 , overdisp ~ Spp -1 ), 
           optimizer = "optimx", method="L-BFGS-B", start = list( Tmax=rep(17,28) , Hmax = rep(0.5,1) , b = rep(0.5,14), overdisp = rep(0.5,14) ) ,
	     lower = c(Tmax=rep(0.1,28) , Hmax = rep(0.5,1) , b = rep(0.01,14), overdisp = rep(0.01,14)),
		 upper =  c(Tmax=rep(200,28), Hmax = rep(10,1), b = rep(3.0,14) , overdisp = rep(10.0,14)),
		 control=list( trace = 3, maxit = 400, factr = 1e12), itnmax =100  )