ans = 
Derive full formula from fitlm(X,y,"quadratic")
Show older comments
I have X which is a 258x6 double and y which is a 258x1 double. These are used to fit a model utilizing fitlm. I found the quadratic to yield the best results.
X = outTableNumeric{:, corr_ids}; % Using the most correlated features
y = outTableNumeric.target; % Target variabel
% Fit a linear regression model
mdl = fitlm(X, y, "quadratic");
However when I try to derive a formula to be deployed outside matlab I find it extremely obscure how I'm to combine the quadratic terms with the coefficients and so on?
I tried to make the formula programmatically but that is not working when I test the result on some row of X.
coefficients = mdl.Coefficients.Estimate;
% Extract the simple notation of the formula
intercept = coefficients(1);
formula = sprintf('%.2f', intercept);
for i = 2:length(coefficients)
if coefficients(i) ~= 0
feature_indices = find(strcmp(mdl.CoefficientNames, mdl.CoefficientNames{i}));
feature_name = mdl.CoefficientNames{feature_indices(1)};
formula = sprintf('%s + %.2f * %s', formula, coefficients(i), feature_name);
end
end
strrep(formula,':','*');
disp('Simple notation of the formula:');
disp(formula);
Can you please help me extract the full equation?
Accepted Answer
John D'Errico
on 23 Jul 2024
Edited: John D'Errico
on 23 Jul 2024
X = rand(258,6);
y = rand(258,1);
mdl = fitlm(X, y, "quadratic")
There are 6 variables, here called x1...x6, corresponding to the 6 columns of your array. And 28 terms in the model. You can see them listed above, as represented by the linear regression model.
If you are not sure how to use the result of fitlm, here are the tools that can be applied:
methods(mdl)
And of course, you can perform a prediction at any point as:
predict(mdl,X(1,:))
If you do extract the coefficients themselves from the model, do NOT extract them only as 5 digit approximations. Do NOT just use the numbers you see under the extimate column. (That is perhaps the most common mistake we see made in MATLAB, thinking that because we see a number written as 0.12097, that is the actual number. WRONG. The actual number is a double precision number. Use it properly.) So we can extract the coefficients as the vector:
coeffs = mdl.Coefficients.Estimate;
format long g
coeffs(1) % the constant term...
13 Comments
Torsten
on 23 Jul 2024
Do you know why linear terms appear for the coefficients, but not in the regression model ?
John D'Errico
on 23 Jul 2024
Edited: John D'Errico
on 23 Jul 2024
Interesting. I did not notice that. I just assumed the linear terms were listed someowhere to the right in that line, but they are not. Looks like a bug in the display utility for a fitlm object. For example, if I do this:
X = rand(258,6);
y = rand(258,1);
mdl = fitlm(X,y)
We see the constant term listed, plus all the linear terms, so 7 terms in total.
mdl2= fitlm(X,y,'quadratic')
If I tell it to use a quadratic model, then it knows there are 28 terms. That corresponds to a constant term, 6 linear terms, 6 squared terms, and 15 crossproduct terms. But the line that says linear regression model seems to have dropped the linear terms in the display.
If I call it with only one independent variable, it shows a correct model:
fitlm(rand(10,1),rand(10,1),'quadratic')
However, when I do the same with two or more independent variables, it seems to be dropping the linear terms from that display.
fitlm(rand(10,2),rand(10,1),'quadratic')
Again, it does know the terms in the model. They can be seen as:
S = struct(mdl2);
S.CoefficientNames
Just that the line in the display seems to be incorrect for problems with 2 or more independent variables. So a bug, but nothing that ever triggered a bug report, because nothing ever actually uses that line from the display. I can post it as a bug to tech support. I would guess the problem lies somewhere in the sub-functions of:
.../MATLAB_R2023a.app/toolbox/stats/classreg/@LinearModel/LinearModel.m
Edit: I've reported this as a bug. However, it is PURELY a bug in the display method for a fitlm object. Since there appear to be no issues in performance due to that bug, that makes it likely they won't put a high priority on fixing it. Of course, since it is surely something simple to fix, I'd not be surprised to see it fixed soon in some new release.
That is a good question from @Torsten, about why the linear terms do not appear in the regression model. It looks to me like a mistake in the console display of the model. The linear terms are being used in the model predictions (see below), even though the console display of the model does not include the linear terms. Simple demonstration of this is below.
y=rand(50,1); X=rand(50,2);
m=fitlm(X,y,'quadratic')
Xa=[ones(50,1),X,X(:,1).*X(:,2),X.^2];
b=m.Coefficients.Estimate;
ypred=Xa*b;
disp([y(1:3),m.Fitted(1:3),ypred(1:3)])
The values m.Fitted and ypred (columns 2 and 3 above) are the same, which shows that the regression model is in fact using the linear terms to compute the Fitted values. Even though the linear regression model displayed on the console, by fitlm(), does not include the linear terms.
John D'Errico
on 23 Jul 2024
Edited: John D'Errico
on 23 Jul 2024
Yes. I have already reported it as a bug to support. However, it is a bug of little real importance, just in that one line of the display.
John D'Errico
on 24 Jul 2024
It seems that behavior is part of Wilkinson notation.
Though I did not see the point clearly made in doc fitlm.
Jacob Emil Elbek
on 24 Jul 2024
You can reconstruct the equation as
fun = @(x) mdl.Coefficients.Estimate(1) + mdl.Coefficients.Estimate(2)*x(1)+mdl.Coefficients.Estimate(3)*x(2) ...
and so on and call the function with a 6x1 input argument.
Jacob Emil Elbek
on 24 Jul 2024
X = rand(258,2);
y = rand(258,1);
mdl = fitlm(X, y, "quadratic")
p = mdl.Coefficients.Estimate
fun = @(x)p(1)+p(2)*x(1)+p(3)*x(2)+p(4)*x(1)*x(2)+p(5)*x(1)^2+p(6)*x(2)^2;
format long
ypred = predict(mdl,[1 6])
fun([1 6])
Unfortunately, they make it more difficult than they should, at least in my opinion, because this was not well documented in the help, AND because they don't have a simple way to extract the explicit order of the terms.
You need to read about
So, when you see this:
Linear regression model:
y ~ 1 + x1*x2 + x1*x3 + x1*x4 + x1*x5 + x1*x6 + x2*x3 + x2*x4 + x2*x5 + x2*x6 + x3*x4 + x3*x5 + x3*x6 + x4*x5 + x4*x6 + x5*x6 + x1^2 + x2^2 + x3^2 + x4^2 + x5^2 + x6^2
the linear terms are all implicitly in the model too. Unfortunately, that does not help you here a lot. Instead, you can find the sequence of terms like this:
X = randn(100,3);y = rand(100,1);
mdl = fitlm(X,y,'quadratic')
mdl.CoefficientNames
As you can see, what they call the linear regression model, is, while valid, based on the Wilkinson notation, it does not lead easily to reconstructing the model with the coefficients in the proper order. However, mdl.CoefficientNames does do that, at least more directly.
@John D'Errico, Thank you for explaining that fitlm() uses Wilkinson notation to display the formula.
If using the quadratic option, then it is not hard to know the right order for the coefficients. As you point out, they are in mdl.CoefficientNames. Or you can simply add the linear terms, in order, after the initial "1". Example.
X=rand(50,3); y=rand(50,1);
mdl=fitlm(X,y,'quadratic');
Display the formula (Wilkinson notation)
mdl.Formula.LinearPredictor
The full formula, with the terms in the right order, is
fullFormula=['1 + x1 + x2 + x3 + ',mdl.Formula.LinearPredictor(5:end)]
If the model order is higher than 2 in some or all predictors, then the simplest approach for the full formula is to actually use mdl.CoefficientNames.
Jacob Emil Elbek
on 24 Jul 2024
John D'Errico
on 25 Jul 2024
Edited: John D'Errico
on 25 Jul 2024
After some discussion with tech support, I learned a resolution, though I did not find it myself. The documentation for LinearModel needs to be enhanced IMHO.
X=rand(50,3); y=rand(50,1);
mdl=fitlm(X,y,'quadratic')
exponents = mdl.Formula.Terms(:,1:end-1)
Note that I needed to drop the last column. But this is essentially a full description of each term in the model. It identifies the exponents of each variable that would be necessary to reconstruct the model. How can we use this? First of all, we should see it does work properly.
syms x1 x2 x3
vpa(mdl.Coefficients.Estimate.'*prod([x1,x2,x3].^exponents,2),4)
For a single point newX, as a row vector, we can just do this:
mdlfun = @(newX) mdl.Coefficients.Estimate.'*prod(newX.^exponents,2);
mdlfun([1 2 3])
And for comparison, we see:
predict(mdl,[1 2 3])
With just a little more effort, we could even vectorize the function.
More Answers (0)
Categories
Find more on Linear Predictive Coding in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
