plot.lrsurf <-
function(Y,X,b0=seq(-2,2,by=0.01),b1=seq(-2,2,by=0.01),ncl=20){
lrs<-outer(b0,b1,"+")*0
for(i in 1:length(Y)){
lpv<-outer(b0,b1*X[i],"+")
lrs<-lrs+lpv*Y[i]-log(1+exp(lpv))
}
image(b0,b1,lrs,xlab=expression(beta[0]),ylab=expression(beta[1]))
contour(b0,b1,lrs,add=T,nlevels=ncl)
}