Nonlinear regression with categorical predictor?

Suppose I have the following equation

 y = (k1|x1)*(k2*x2)*(k3*x3^(1/k4))

where k1 is a coefficient depending on variable x1. Both k2 and k3 are also coefficients, and k4 is a power term, which all the ks are not known to me.

x1 is categorical, say adult male, young male, adult female and young female, and x2 and x3 are numeric, say x2 = weight 75kg, 62kg, 89kg... and x3 = height 180 cm, 172 cm, 170 cm...

Anyone knows how to perform a regression for such a combination of data to find all the ks? and eventually the model has two values for k1, for example: if x1 = male, k1 = 2.5; if x1 = female, k1 = 1.5.

[the cyclist]

I haven't thought about how to model this whole thing, but the term

k1*x1

is problematic, I think, when x1 is categorical. For example, what does "6 times male" mean?

Since you didn't mention explicitly, I assume that x2 and x3 are interval data?

[wesleynotwise]

Ah. Good question. Let me modify my the first term and my question. x2 and x3 are interval data, for example, weight and height in this case.

Answer 1

[Michelangelo Ricciulli]

Ok, a very simple way to do it is the following.

You can use the function fminunc, that finds the minimum value of something. What is this something? You want your model to predict very well the y value, so the mean square error between y and the model is what you want to minimize. Let's define this function as errFunc depending on a vector param:

errFunc=@(param) mean((y - k1(x1)*(param(1)*x2)*(param(2)*x3.^(1/param(3)))).^2);

This works if you already have in your environment the data x1, x2,x3,y and also the function k1 that returns the right coefficient based on k1.

Then, you just need to call fminunc with the function you just created and a guess of the 3 values you are searching (let's just put random numbers as guess)

   fminunc(errFunc,randn(3,1))

this will output the value of param you are searching for.

Probably you'll need to add another another coefficient, let's call it param(4), that is summed to your model to better fit the data.

[wesleynotwise]

Hi. At the moment, I have all the variable data, but all the k terms in my example (or parameters in your example) are not known to me (that's why regression/ model is needed).

Regarding to the 4th parameter, do you mean the k0 in the following equation

y = k0 + (k1|x1)*(k2.x2)*(k3.x3^(1/k4))

[Michelangelo Ricciulli]

Yes, I meant exactly that. So also the k1 is not known to you? Then I would do it like this: first I'll transform my x1 in a binary vector.

 x1new=(x1=='Male');

Then, you can write the first part of the model as depending on 2 new variables k1 and k1_new like this

 ...(k1_male*x1_new + k1_female*(1-x1_new))...

This should be equivalent to having a constant k1 that changes its value according to the sex. In fact when x_new will be 1, the constant for female will be multiplied by (1-x1)=0.

[wesleynotwise]

Interesting.... Two questions: (1) Is it possible to drop the k0 from the regression? (2) What if I have 4 categorical arrays? eg: adult male, young male, adult female and young female?

[Jan]

@weskeynotwise: You have flagged your own comment as spam. I think the question does not conflict with the terms of usage and removed the flag. Well, I think you had another intention.

[wesleynotwise]

@Jan Simon: When I posted my comment, the system recognised it as a 'potential spam', and asked me to either flag it for review or contact the administration. I chose the former.

[Michelangelo Ricciulli]

(1) I think it is possible, but I'm not sure whether the fit would be as good. Anyway it depends on your data, you can give it a try.

(2) If you have 4 categories, I would transform the vector x1 in 3 (not 4) binary different vector. So a vector x1_ym that contains 1 only if x1 has 'male young', a vector x1_ya with 1 only if x1 has 'male adult', and a vector x1_wy with 1 if x1 has 'female young'. Then I would add 4 parameters like this:

...( param(1)*x1_ym+param(2)*x1_am+param(3)*x1_yw+param(4))...

By the way, fminunc returns two things, a vector containing the variable, and the corresponding value of the error. When you try different approaches you can check the value of the error: the lower, the better.

[wesleynotwise]

Thanks for your quick reply. (1) I think in my case k0 has less meaning, because if all the variables = 0 (no input), the output of y should be 0, instead of a fixed constant of k0. (2) This is getting slightly complicated now. Why the vector is 3 but not 4? In the example, 3 vectors are defined for adult male, young male and young female, but what about the adult female?

Thanks in advance!

[Michelangelo Ricciulli]

Usually, with N categorical variables, you need to use N-1 binary variables.

The reason is simple: if you know that someone is not a young-male, nor a adult-male, nor a young female, it will certainly be an adult female. So, it's somehow already included in your data, and including it twice could make things more difficult for the optimization algorithm.

[wesleynotwise]

Cool :) Let me recap. And the code should look like this?

tbl = table(x1, x2, x3, y);
fit = @(k,x) (k1*x1_am+k2*x1_ym+k3*x1_af)*(k4.x2)*(k5.x3^(1/k6));
k0 = [1 1 1 1 1 1]; % just assign any random number
mdl = fitnlm(tbl,fit,k0)

[Michelangelo Ricciulli]

There are some adjustments to be done. (1) table should already contain the modified data. So you need to first create x1_am, x1_ym and x1_af. Then issue:

 tbl = table(x1_am, x1_ym, x1_af, x2, x3, y);

(2) k and x in your fit function are respectively a vector and a matrix. So you need to define your fucntion depending on these variables:

 fit = @(k,x) (k(1)*x(:,1)+k(2)*x(:,2)+k(3)*x(:,3))*(k(4)*x(:,4))*(k(5).x(:,5).^(1/k(6)));

(3) you will need a k(7) in the first part, otherwise it will always be 0 when you have a female adult. Like this

fit = @(k,x) (k(1)*x(:,1)+k(2)*x(:,2)+k(3)*x(:,3) +k(7))*(k(4)*x(:,4))*(k(5).x(:,5).^(1/k(6)));

so you'll also need to add another value to the k0, since now there are seven parameters to optimize.

[Michelangelo Ricciulli]

Check the other answer, it makes a very good point. Sorry if I didn't notice that before

[wesleynotwise]

I tried your suggestion for a small sample size, and... IT WORKS!!! You've brighten up my day.... Well, I receieved the following error message though, I assume that was due to my data size, or the variables I used, which should be easy to resolve.

Warning: The Jacobian at the solution is ill-conditioned, and some model parameters may not be estimated well
(they are not identifiable).  Use caution in making predictions. 
> In nlinfit (line 376)
In NonLinearModel/fitter (line 1123)
In classreg.regr.FitObject/doFit (line 94)
In NonLinearModel.fit (line 1430)
In fitnlm (line 94)

Answer 2

[Ilya]

Unless I misunderstood your dot notation, the problem is ill-defined. It has an infinite number of solutions. Rewrite it in this form:

y = ((k1|x1)*k2*k3) * x2*x3^(1/k4)

Observe that you can only fit for ((k1|x1)*k2*k3), but not for separate coefficients k1, k2 and k3.

Generally, the best way to handle multiplicative models is to turn them into additive models by taking the log of both sides. If you do that, you get

q = (c1|x1) + c4*z3

where

q = log(y) - log(x2)

c1 = log(k123|x1)

k123|x1 = (k1|x1)*k2*k3

c4 = 1/k4

z3 = log(x3)

This model can be easily fitted by fitlme. (If you have only two levels in x1, you can easily get away without fitlme.) The formula would be something like 'q ~ (1|x1) + z3'.

[wesleynotwise]

Hi IIya. Sorry I didn't see your message early. Yes. The dots and asterisks mean multiply in my example. Just realise that mistake. I'll correct them. The expression should look like this

   y = (k1|x1)*(k2*x2)*(k3*x3^(1/k4))

where k1 depends on x1, which is a categorical variable, eg: adult male, young male, adult female and young female. So, the equation will have 4 different k1. I think solving my problem using linear regression is another possible method, which I want to give it a try! So, If converting the equation to log term, it becomes

q = (c1|x1) + c4*z3

where q = log (y) - log (x2);
c1|x1 = log (k1|x1) + log (k2) +log (k3);
c4 = 1/k4;
z3 = log (x3);

after running the linear regression, the matlab is likely to give me values for c1 for different x1 and c4. Does that mean I still have to do a separate regression to find k1, k2 and k3? My intention is to have the equation in the original form (non-log), so that one just need to feed in x1, x2 and x3 to find y in the future. Any thoughts?

[Ilya]

You can't separate k1, k2 and k3. The form of your problem does not allow that. You can feed your problem into some minimizer and get a numeric answer (likely with some warnings about ill-conditioning). But the numeric answer such as {k1a,k1b,k2,k3} (a and b are used to enumerate categorical levels) is not in any way different from {k1a,k1b,k2/2,2*k3} or {2*k1a,2*k1b,k2/2,k3} etc. If you run fitlme, you will get K values for the intercept c1=log(k1*k2*k3), where K is the number of categorical levels in x1, plus one linear coefficient for c4. That's all you can do with this problem.

[wesleynotwise]

Ah. If that's the case, the equation can at most be re-written to this

log (y) = (c1|x1) + c4*log(x3) + log (x2)

This means once one feeds all the variable data of x1, x2 and x3, one needs to inverse the log of y in order to obtain real y value?