Compare MLR model to model $Y_i = (\beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \epsilon_i)^{\beta_3}$?

Question

If I have theoretical reasons to suppose the data might be fit with an unusual equation such as the following:

$$Y_i = (\beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \epsilon_i)^{\beta_3}$$

Can I use Ordinary Least Squares Multiple Linear Regression after a transformation to estimate parameters $\beta{_0,_1,_2,_3}$? If yes, what transformation?

If not, is there some specialized package in R (and brief reading) that might help me compare the fit and residuals from this model against a more typical MLR model?

Thanks.

Example Code:

## while I can run "nls," I cannot get $\epsilon$ inside parentheses nor
## can I have four BETAs

var1 <- rnorm(50, 100, 1)
var2 <- rnorm(50, 120, 2)
var3 <- rnorm(50, 500, 5)

## make a model without $\beta_1$ and $\beta_2$ and with $\epsilon_i$ on outside
nls(var3 ~ (a + var1 + var2)^b, start = list(a = 0.12345, b = 0.54321))

Nonlinear regression model
  model: var3 ~ (a + var1 + var2)^b
  data: parent.frame()
   a        b 
 475.5234   0.9497 
 residual sum-of-squares: 1365

Number of iterations to convergence: 6 
Achieved convergence tolerance: 8.332e-08

## FAILS with exponent on left-hand side and $\epsilon$ inside parentheses
nls(var3^(1/b) ~ (a + var1 + var2), start = list(a = 0.12345, b = 0.54321))
Error in eval(expr, envir, enclos) : object 'b' not found

## FAILS with all BETAs
nls(var3 ~ (a + b*var1 + c*var2)^d, start = list(a = 4, b = 1, c = 1, d = 1))
Error in numericDeriv(form[[3L]], names(ind), env) : 
Missing value or an infinity produced when evaluating the model

Answer 1

No (at least not with nls)

From its documentation, nls fits functions of the form $Y_i| \theta, X_i = f(\theta, X_i) + \epsilon$ (and is the MLE in the case that $\epsilon$ is iid Normal), so your relationship is not in the non-linear least squares class.

Let's see if we can describe the distribution $Y$ might follow. Let $Z_i = \beta_0+\beta_1 x_{1i} + \beta_2 x_{2i} + \epsilon_i$ Given that $\epsilon_i$ is $N(0, 1)$, then $Z_i \sim N(\beta_0+\beta_1 x_{1i} + \beta_2 x_{2i}, 1)$. If $\beta_3 = 2$ then for example, we could have that $Y_i$ is non-central $\chi^2_1$.

Yes (using Box-cox transformations)

If $Y_i = Z_{i}^{\beta_3}$ is a one-to-one transformation (ie, at a minimum, $\beta_3$ is not even) then you have just rediscovered the box-cox family of transformations: $$ Y(\lambda) = \begin{cases} (\lambda Z + 1)^{1/\lambda}, \lambda >0 \\ e^Z, \lambda = 0 \end{cases}, $$ which clearly includes the scenario you describe. Classically, $\lambda$ is estimated through the profile likelihood, ie, plugging in different values of $\lambda$ and checking the RSS to the least-squares fit. An Analysis of Transformations Revisited (1981) appears to give a good review of the theory. The function boxcox in the package MASS does such an estimation. If $\beta_3$ is a parameter of interest rather than a nuisance you may need to do something more sophisticated.

Answer 2

I think Andrew M have given a good answer; I just want to make a few related points.

As Andrew M indicates you can't do the model as is directly with nonlinear least squares; however, you can fit this closely related model with nonlinear LS:

$Y_i = (\beta_0 + \beta_1x_{1i} + \beta_2x_{2i})^{\beta_3} + \epsilon_i$

This might not seem much use, but it would have value in obtaining an initial estimate of $\beta_3$ to get a good starting point for optimization of the actual model (whether performed directly, or via Box-Cox).

Note also that if $Y$ is strictly positive, you can consider this transformation:

$\log(Y_i) = \beta_3 \log(\beta_0 + \beta_1x_{1i} + \beta_2x_{2i} + \epsilon_i)$

Again, a slight modification (pulling the error term outside the parentheses) allows nonlinear least squares fitting. You could then reweight using the resulting estimate of $\beta_3$ to improve the estimates. The only difficulty would be if you hit a situation where the fitted value inside the log wasn't strictly positive.

[If you're prepare to consider Weibull regression (that is, where the Y's are Weibull with mean dependent on the X's), you might find that you can do something useful with that. It would change the form of the relationship with the x's however. A related approach would be that given a value for $\beta_3$ you could consider transforming $Y$ ($Y^*=Y^{1/\beta_3}$)and fit an exponential GLM with identity link to $Y^*$ rather than a Gaussian. This would again correspond to a Weibull model for $Y$, but with the parameters entering in the way you suggest). This could be done over a grid of $\beta_3$ values to maximize the likelihood for it.]