[R] constrained logistic regression: Error in optim() with method = "L-BFGS-B"

[ashky]

Dear R Users/Experts,

I am using a function called logitreg() originally described in MASS (the
book 4th Ed.) by Venebles & Ripley, p445. I used the code as provided but
made couple of changes to run a 'constrained' logistic regression, I set the
method = "L-BFGS-B", set lower/upper values for the variables.

Here is the function,

	logitregVR <- function(x, y, wt = rep(1, length(y)), intercept = T, start =
rep(0.1, p), ...)
	{

	#-------------------------------------------#
	# A function to be minimized (or maximized) #
	#-------------------------------------------#
 	 fmin <- function(beta, X, y, w) 
		{
			p <- plogis(X %*% beta)
			-sum(2*w*ifelse(y, log(p), log(1-p)))		# Residual Deviance: D =
-2[Likelihood Ratio]
			#print(-sum(2*w*ifelse(y, log(p), log(1-p))))
		}

	#----------------------------------------------------------------------#
	# A function to return the gradient for the "BFGS", "L-BFGS-B" methods #
	#----------------------------------------------------------------------#
 	 gmin <- function(beta, X, y, w) 
		{
			eta <- X %*% beta; p <- plogis(eta)
			-2*matrix(w*dlogis(eta)*ifelse(y, 1/p, -1/(1-p)), 1) %*% X		# Gradient
			#print(-2*matrix(w*dlogis(eta)*ifelse(y, 1/p, -1/(1-p)), 1) %*% X)
  		}

  	 if(is.null(dim(x))) dim(x) <- c(length(x), 1)
  	 dn <- dimnames(x)[[2]]
  	 if(!length(dn)) dn <- paste("Slope", 1:ncol(x), sep="")
  	 p <- ncol(x) + intercept # 1 + 1
  	 if(intercept) {x <- cbind(1, x); dn <- c("(Intercept)", dn)}
  	 if(is.factor(y)) y <- (unclass(y) != 1)

  	 fit <- optim(start, fmin, gmin, X = x, y = y, w = wt, method =
"L-BFGS-B", lower = c(-Inf, 0.05), upper = c(Inf, Inf), ...)
  	 #fit <- optim(start, fmin, gmin, X = x, y = y, w = wt, control =
list(trace = 6), method = "L-BFGS-B", lower = c(-Inf, 0.05), upper = c(Inf,
Inf), ...)

  	 names(fit$par) <- dn
  	 cat("\nCoefficients:\n"); print(fit$par)
  	 cat("\nResidual Deviance:", format(fit$value), "\n")
  	 if(fit$convergence > 0) cat("\nConvergence code:", fit$convergence,
"\n")
	 return(list(fitpar = fit$par))
  	 invisible(fit)
	}

And here is my data,

# ----------------------
# Dose     0 1 emp.prop 
# ----------------------  
# 1           3 0  0.000
# 2           1 0  0.000
# 3.3        4 0  0.000
# 5.016     3 0  0.000
# 7.0224    4 0  0.000
# 9.3398    2 0  0.000
# 12.4219  2 0  0.000
# 16.5211 47 1  0.021
# 21.9731  1 1  0.500
# ---------------------- 

Clearly, I can use the glm() function with family = binomial (or) the lrm()
function. But my interest is to constraint the slope parameter (beta). If
you see the above above by Venebles and Ripley, I chose method = "L-BFGS-B"
and set lower and upper values for my two parameters. Here is the code and
the error I am getting, 

y <- c(rep(0,3), rep(0,1), rep(0,4), rep(0,3), rep(0,4), rep(0,2), rep(0,2),
rep(0,47), rep(1,1), rep(0,1), rep(1,1))
x <- c(rep(1,3), rep(2,1), rep(3.3,4), rep(5.016,3), rep(7.0224,4),
rep(9.3398,2), rep(12.4219,2), rep(16.5211,48), rep(21.9731,2))

length(x)
length(y)

	#------------------#
	# Method 1 - Works #
	#------------------# 
	glm(y~x, family = binomial)$coefficients

	# summary(glm(y~x, family = binomial)) 

	#------------------#
	# Method 2 - Works #
	#------------------# 
             library(Design)
	lrm(y ~ x)$coef	

	#-------------------#
	# Method 3 - Works #
	#-------------------# 
	lgtreg <- function(beta, x, y)
	{ 
		eta <- exp(beta[1] + beta[2]*x)
		p <- eta/(1+eta)
		return(-sum(ifelse(y,log(p),log(1-p)))) 	# This is the log-likelihood
	}
           optim(c(0,0), lgtreg, x = x, y = y, method = "BFGS", hessian =
TRUE)$par  

	#---------------------------------------------#
	# Method 4 - Veneble and Ripley's - Gives error #
	#--------------------------------------------#
	logitregVR(x, y)

I am getting this error,

Error in optim(start, fmin, gmin, X = x, y = y, w = wt, method = "L-BFGS-B", 
: L-BFGS-B needs finite values of 'fn'

I am really not sure whats causing this error, can someone please help me.
Any comments/replies highly appreciated,

Ashky
-- 
View this message in context: http://www.nabble.com/constrained-logistic-regression%3A-Error-in-optim%28%29-with-method-%3D-%22L-BFGS-B%22-tp19709337p19709337.html
Sent from the R help mailing list archive at Nabble.com.