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

Question

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

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

Hay una "rápida" manera de hacer esto con pre-existente de las funciones R o hay un nombre para un modelo como este? Me doy cuenta de que se puede modificar el método de Newton-Raphson algoritmo utilizado para la regresión logística, pero esto es un montón de teórico y el trabajo de codificación y estoy buscando un atajo.

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

Answer 1

Esto es bastante fácil de hacer con la optim función en R. Mi entendimiento es que desea ejecutar una regresión logística, donde y es el binario. Simplemente escribe la función y, a continuación, pegarla en optim. A continuación está el código que no se ha ejecutado (pseudo có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

Aviso que tu.la diversión es el negativo de un registro de probabilidad de la función. Así optim es maximizar el registro de la probabilidad (por defecto optim minimiza todo lo que es la razón por la que hice la función negativa). Y si no es binario ir a http://fisher.osu.edu/~schroeder.9/AMIS900/ch5.pdf para multinomial condicional y la función de las formas en los modelos logit.

Answer 2

La respuesta anterior es correcta. Para referencia, aquí hay algunos elaborados de trabajo código R para calcular. Tengo que tomar la libertad de agregar una intercepción, porque probablemente lo quiero uno de esos.

## 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, la construcción de un registro de probabilidad de la función a maximizar, por medio de la dbinom porque es allí, y sumar 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 máxima verosimilitud. Yo no he molestado a ofrecer un gradiente o elegir un método de optimización, pero es posible que desee hacer tanto.

## 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 mira los resultados. El ML estimaciones de los parámetros y asintótica SEs son:

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

que debe ser

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

o hay un bug (que siempre es posible).

Las habituales advertencias acerca de Hesse-derivados de los errores estándar se aplican.