8
$\begingroup$

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.

$\endgroup$
12
  • 1
    $\begingroup$ Can you explain the situation in which you'd want to do that ever? I'm just curious. I personally haven't seen this and would probably have to code it up using constrained optimization. $\endgroup$ – wolfsatthedoor Dec 4 '14 at 18:12
  • 1
    $\begingroup$ I cant get into the specifics but the reason I would want to do this is basically for statistical power. If I believe this relationship exists then forcing the model to have this relationship will get me closer to the true $\beta$ value. As for getting point estimates on $\beta$ , thats pretty easy using optim() or some other optimizer in R to maximize the likelihood. But I need standard errors on these guys. $\endgroup$ – TrynnaDoStat Dec 4 '14 at 18:18
  • 3
    $\begingroup$ This one is not as hard as it might look: it's pretty easy to obtain the MLE from first principles, because all you have is a single parameter. Write down the log likelihood and differentiate it with respect to $\beta$. Find the zero(s). That would be done numerically. If you have any trouble starting the search, fit the usual two-parameter model and use (say) the coefficient of $x_2$ as the starting value. $\endgroup$ – whuber Dec 4 '14 at 18:56
  • 2
    $\begingroup$ The NLIN procedure in SAS can fit a nonlinear regression formula such as this and will output coefficient stnd errors. Are you married to R or is there some flexibility? $\endgroup$ – RobertF Dec 4 '14 at 19:10
  • 1
    $\begingroup$ Getting the SEs isn't any harder: this is the easy case of the standard theory, where there is only one parameter. But given the nonlinear nature of its dependence on the parameters, I would be reluctant to relie either on the approximate theory or on brute-force numerical optimization: plot the likelihood function itself, at least in a few cases, so you can understand its qualitative behavior near the peak. $\endgroup$ – whuber Dec 4 '14 at 22:58
13
$\begingroup$

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.

$\endgroup$
4
  • 1
    $\begingroup$ Appreciate the response Zachary! However, this won't work as I need standard errors on my estimates. I'm thinking of combining bootstrapping and optim() but would really prefer a non-bootstrapping method if possible. Modifying Newton-Raphson would be much more satisfying but much harder to implement. $\endgroup$ – TrynnaDoStat Dec 4 '14 at 22:04
  • 3
    $\begingroup$ Maybe I don't understand, the standard error of the estimate comes from a function of the maximum likelihood hessian evaluated at the estimates. The way you wrote your function you only have one parameter. You can certainly bootstrap the above code to obtain standard errors too. $\endgroup$ – Zachary Blumenfeld Dec 4 '14 at 22:19
  • 1
    $\begingroup$ @ZacharyBlumenfeld I understand what you're saying now! I was confused because my understanding of asymptotic theory of MLE was that our observations must be iid (the mean certainly varies here so my observations are not iid). However, I see now that observations do not have to be iid under certain regularity conditions (en.wikipedia.org/wiki/Maximum_likelihood#Asymptotic_normality). Now I just need to double check my situation meets the regularity conditions. Thanks again! $\endgroup$ – TrynnaDoStat Dec 5 '14 at 15:09
  • 1
    $\begingroup$ Note: If $Y = X\beta + \epsilon$ then it is $\epsilon$ that might be assumed iid, not necessarily $Y$. $\endgroup$ – conjugateprior Dec 5 '14 at 17:46
2
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.