Maximize nonlinear regression function in R

Question

Given a linear model obtained from the function call reg = lm(...), how can you find the coefficients that maximize the obtained regression function?

I'm aware of the function optim(...), but it requires a function as an input. I haven't figured out how to extract this from the regression model.

It should be noted that I'm using non-linear terms in my regression analysis (squared variables, to be precise).

In other words, by regression function looks like

y_hat = kx_11*x_1+kx_12*x_1^2 + kx_21*x_2+kx_22*x_2^2 + ...

What do you mean by "coefficients that maximize the obtained regression function"? — thc
I just updated the post for clarity. What I'm trying to achieve corresponds to taking the partial derivatives with respect to all variables and solving the set of equations in which all partial derivatives equal zero. — Eivind
You didn't clarify at all. Taking your question literally, the regression function will not be maximized because you can always increase the magnitude of your coefficients without bound. I.e., "setting your coefficients to infinity" will produce a "maximum." I assume that's not what you want to do, so you need to clarify more. — thc
No, because the regression has quadratic terms, as I pointed out. x - x^2, for instance, does not have a maximum at infinity. — Eivind
Finding the maximum of a formula is a different question than finding the regression of a model (i.e., the goodness of fit). Input into the lm function would optimize the fit of y=Ax - Bx^2. In your case, you have no y values, so it doesn't make sense to use an R formula. — thc

danielson danielson · Accepted Answer · 2017-01-04T02:18:51

Here is a quick example to demonstrate 1 way. Use predict() on the lm object to create your function. fxn() is a little messy since I don't have your exact data, but you should get the idea.

#set up dummy data
x1 = -10:10
x2 = runif(21)
y = -x1^2 + x1 - 10*x2^2 + runif(21)*.1 
data = data.frame(y= y, x1=x1, x2=x2)

#fit model
m = lm(data=data, y ~ x1 + I(x1^2) + I(x2^2))

#define function that returns predicted value
fxn = function(z){
    z = as.data.frame( t(z) )
    colnames(z) = colnames(data)[-1]
    predict(m, newdata=z)
}

optim(c(0,0), fxn, control=list(fnscale=-1)) #maximizes fxn

$par
[1]  4.991601e-01 -3.337561e-06

$value
[1] 0.3153461

$counts
function gradient 
      65       NA 

$convergence
[1] 0

$message
NULL

Maximize nonlinear regression function in R

1 Answers