Existe-t-il un moyen de forcer une relation entre les coefficients de régression logistique?

Je voudrais spécifier un modèle de régression logistique où j'ai la relation suivante:

$E[Y_i|X_i] = f(\beta x_{i1} + \beta^2x_{i2})$ où est la fonction logit inverse.$f$

Existe-t-il un moyen "rapide" de le faire avec des fonctions R préexistantes ou existe-t-il un nom pour un modèle comme celui-ci? Je me rends compte que je peux modifier l'algorithme de Newton-Raphson utilisé pour la régression logistique, mais c'est beaucoup de travail théorique et de codage et je cherche un raccourci.

EDIT: obtenir des estimations ponctuelles pour est assez facile en utilisant optim () ou un autre optimiseur dans R pour maximiser la probabilité. Mais j'ai besoin d'erreurs standard sur ces gars-là.$\beta$

C'est assez facile à faire avec la fonction optim dans R. Je crois comprendre que vous voulez exécuter une régression logistique où y est binaire. Vous écrivez simplement la fonction, puis la collez dans optim. Vous trouverez ci-dessous un code que je n'ai pas exécuté (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

Notez que your.fun est le négatif d'une fonction de vraisemblance de journal. Donc, optim maximise la vraisemblance du journal (par défaut, optim minimise tout ce qui explique pourquoi j'ai rendu la fonction négative). Si Y n'est pas binaire, allez à http://fisher.osu.edu/~schroeder.9/AMIS900/ch5.pdf pour les formes de fonctions multinomiales et conditionnelles dans les modèles logit.

La réponse ci-dessus est correcte. Pour référence, voici un code R élaboré élaboré pour le calculer. J'ai la liberté d'ajouter une interception, car vous en voulez probablement une.

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

Construisez maintenant une fonction de vraisemblance logarithmique à maximiser, en l'utilisant ici dbinomparce qu'elle est là, et en additionnant les résultats

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

et ajuster le modèle par maximum de vraisemblance. Je n'ai pas pris la peine d'offrir un dégradé ou de choisir une méthode d'optimisation, mais vous voudrez peut-être faire les deux.

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

Jetez maintenant un œil aux résultats. Les estimations des paramètres ML et les SE asymptotiques sont:

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

qui devrait être

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

ou il y a un bug (ce qui est toujours possible).

Les mises en garde habituelles concernant les erreurs standard dérivées de la Hesse s'appliquent.