Is there a way to force a relationship between coefficients in logistic regression?

Question

I would like to specify a logistic regression model where I have the following relationship:

$E[Y_i|X_i] = f(\beta x_{i1} + \beta^2x_{i2})$ where $f$ is the inverse logit function.

Is there a "quick" way to do this with pre-existing R functions or is there a name for a model like this? I realize I can modify the Newton-Raphson algorithm used for logistic regression but this is a lot of theoretic and coding work and I'm looking for a short cut.

EDIT: getting point estimates for $\beta$ is pretty easy using optim() or some other optimizer in R to maximize the likelihood. But I need standard errors on these guys.

Answer 1

This is fairly easy to do with the optim function in R. My understanding is that you want to run a logistic regression where y is binary. You simply write the function and then stick it into optim. Below is some code I didn't run (pseudo code).

#d is your data frame and y is normalized to 0,1
your.fun=function(b)
{

    EXP=exp(d$x1*b +d$x2*b^2)

    VALS=( EXP/(1+EXP) )^(d$y)*( 1/(1+EXP) )^(1-d$y) 
    return(-sum(log(VALS)))
}

result=optim(0,your.fun,method="BFGS",hessian=TRUE)
# estimates 
 result$par
    #standard errors
    sqrt(diag(inv(result$hessian)))
# maximum log likelihood
-result$value

Notice that your.fun is the negative of a log likelihood function. So optim is maximizing the log likelihood (by default optim minimizes everything which is why I made the function negative). If Y is not binary go to http://fisher.osu.edu/~schroeder.9/AMIS900/ch5.pdf for multinomial and conditional function forms in logit models.

Answer 2

The answer above is correct. For reference, here's some elaborated working R code to compute it. I have take the liberty of adding an intercept, because you probably do want one of those.

## make some data
set.seed(1234)
N <- 2000
x1 <- rnorm(N)
x2 <- rnorm(N)
## create linear predictor
lpred <- 0.5 + 0.5 * x1 + 0.25 * x2
## apply inverse link function
ey <- 1/(1 + exp(-lpred))
## sample some dependent variable
y <- rbinom(N, prob=ey, size=rep(1,N))

dat <- matrix(c(x1, x2, y), nrow=N, ncol=3)
colnames(dat) <- c('x1', 'x2', 'y')

Now construct a log likelihood function to maximise, here using dbinom because it's there, and summing the results

## the log likelihood function
log.like <- function(beta, dat){
  lpred <- beta[1] + dat[,'x1'] * beta[2] + dat[,'x2'] * beta[2]**2
  ey <- 1/(1 + exp(-lpred))
  sum(dbinom(dat[,'y'], prob=ey, size=rep(1,nrow(dat)), log=TRUE))
}

and fit the model by maximum likelihood. I haven't bothered to offer a gradient or choose an optimisation method, but you might want to do both.

## fit
res <- optim(par=c(1,1), ## starting values 
             fn=log.like,
             control=list(fnscale=-1), ## maximise not minimise
             hessian=TRUE, ## for SEs
             dat=dat)

Now have a look at the results. The ML parameter estimates and asymptotic SEs are:

## results
data.frame(coef=res$par,
           SE=sqrt(diag(solve(-res$hessian))))

which should be

##        coef         SE
## 1 0.4731680 0.04828779
## 2 0.5799311 0.03363505

or there's a bug (which is always possible).

The usual caveats about Hessian-derived standard errors apply.