¿Hay alguna manera de forzar una relación entre coeficientes en la regresión logística?

Me gustaría especificar un modelo de regresión logística donde tenga la siguiente relación:

$E[Y_i|X_i] = f(\beta x_{i1} + \beta^2x_{i2})$ donde es la función de logit inversa.$f$

¿Existe una forma "rápida" de hacer esto con funciones R preexistentes o hay un nombre para un modelo como este? Me doy cuenta de que puedo modificar el algoritmo de Newton-Raphson utilizado para la regresión logística, pero esto es mucho trabajo teórico y de codificación y estoy buscando un atajo.

EDITAR: obtener estimaciones puntuales para es bastante fácil usando optim () o algún otro optimizador en R para maximizar la probabilidad. Pero necesito errores estándar en estos tipos.$\beta$

Esto es bastante fácil de hacer con la función de optimización en R. Entiendo que desea ejecutar una regresión logística donde y es binario. Simplemente escribe la función y luego la pega en optim. A continuación se muestra un código que no ejecuté (pseudocódigo).

#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

Tenga en cuenta que your.fun es el negativo de una función de probabilidad de registro. Entonces optim está maximizando la probabilidad de registro (por defecto, optim minimiza todo, por eso hice que la función sea negativa). Si Y no es binario, vaya a http://fisher.osu.edu/~schroeder.9/AMIS900/ch5.pdf para obtener formas de funciones condicionales y multinomiales en modelos logit.

La respuesta anterior es correcta. Como referencia, aquí hay un código R de trabajo elaborado para calcularlo. Me he tomado la libertad de agregar una intercepción, porque probablemente quieras una de esas.

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

Ahora construya una función de probabilidad de registro para maximizar, usando aquí dbinomporque está allí y sumando los resultados

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

y ajustar el modelo por la máxima probabilidad. No me he molestado en ofrecer un gradiente o elegir un método de optimización, pero es posible que desee hacer ambas cosas.

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

Ahora eche un vistazo a los resultados. Las estimaciones del parámetro ML y los SE asintóticos son:

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

cual debería ser

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

o hay un error (que siempre es posible).

Se aplican las advertencias habituales sobre los errores estándar derivados de Hesse.