fitlm returns pvalues equal to NaN without zscoring
Show older comments
I can not understand which is the reason why the fitlm using variables without zscoring returns pvalues equal to NaN, whereas this does not happen using zscored variables. Below you can find the code I used:
medians_var1_scored = nanzscore(medians_var1);
medians_var2_scored = nanzscore(medians_var2);
medians_var1_scored_tr = medians_var1_scored';
medians_var2_scored_tr = medians_var2_scored';
tbl=table(medians_var1_scored_tr,medians_var2_scored_tr,'VariableNames', ...
{'var1','var2'});
%build your model
mdl=fitlm(tbl,'var1 ~ var2','RobustOpts','on')
which gives this result:
mdl =
Linear regression model (robust fit):
var1 ~ 1 + var2
Estimated Coefficients:
Estimate SE tStat pValue
_________ ________ ________ _______
(Intercept) 0.028674 0.094595 0.30312 0.76234
var2 -0.072919 0.094998 -0.76758 0.44429
Number of observations: 118, Error degrees of freedom: 116
Root Mean Squared Error: 1.03
R-squared: 0.00584, Adjusted R-Squared: -0.00273
F-statistic vs. constant model: 0.681, p-value = 0.411
If instead I use the original variables, I do:
tbl=table(medians_var1',medians_var2','VariableNames', ...
{'var1','var2'});
%build your model
mdl=fitlm(tbl,'var1 ~ var2','RobustOpts','on')
and obtain:
mdl =
Linear regression model (robust fit):
var1 ~ 1 + var2
Estimated Coefficients:
Estimate SE tStat pValue
__________ __________ ______ __________
(Intercept) 0 0 NaN NaN
var2 8.6386e-13 2.1087e-14 40.966 7.8796e-71
Number of observations: 118, Error degrees of freedom: 117
Root Mean Squared Error: 31.4
R-squared: 0.263, Adjusted R-Squared: 0.263
F-statistic vs. constant model: Inf, p-value = NaN
I can't understand why this happens. You can find attached both the original var1 and var2 variables and the zscored ones. But If I plot both of them, I obtain the same plot (just rescaled). Hence, there shouldn't be a problem in the nanzscore function which is "hand written".
1 Comment
Torsten
on 31 Dec 2025
You forgot to include "nanzscore".
Accepted Answer
d=dir('median*.mat'); % load the files w/o having to type a zillion characters...
for i=1:numel(d)
load(d(i).name)
disp(d(i).name) % figure out which is var1, var2??? show which file
whos('-file',d(i).name) % contains which variable...
end
whos
es=medians_energy_scored_tr; ds=medians_dwi_scored_tr; % shorten names
e=medians_energy_vec.'; d=medians_micro_parameter_vec.';
[all(isfinite(es)), all(isfinite(ds))] % is there a NaN lurking, maybe?
[all(isfinite(e)), all(isfinite(d))]
[min(es) max(es) min(ds) max(ds)] % see what the magnitude of variables are...
[min(d) max(d)], [min(e) max(e)]
Oh....there's likely the problem; d (medians_micro_parameter_vec) is huge in magnitude so the unscaled calculation overflows when computing the sum terms in X'X. The scaling brings everything down to roughly +/-3 range if is normal scaling.
Always examine data before blindly fitting...
plot(ds,es,'x')
figure
plot(d,e,'x')
These appear to have just been randomly generated values---why such a large magnitude was used is anybody's guess here...but there's the problem.
We'll just fit the same data, still large but 10E4 instead of 10E14...
fitlm(d/1E10,e,'RobustOpts','on')
The problem is simply one of the magnitude of the independent variable is so large the sums during computation lose precision leading to an apparently singular design matrix; as the above shows, simply a reduction in the absolute magnitude of the x variable is sufficient; it isn't necessarily there being "magic" in z-scoring.
First we checked to be sure there weren't any NaN that were causing the issue that the function nanzscore() was silently taking care of for us...
If we then adjust by the mean and scale, but only by a power of 10 instead of std(), we get roughly the same result of about a zero intercept, but note the slope is the same only scaled by the additional 10E4 we introduced, but the pValue is also identical.
fitlm((d-mean(d))/1E14,e,'RobustOpts','on')
And just for good measure, we'll use standardized variables...
fitlm(zscore(d),zscore(e),'RobustOpts','on')
This effect of keeping intermediate calculations of the design matrix roughly on the same order as the constant column of ones is one of the prime uses of z-scaling in regression besides putting variables of differing magnitudes on the same scale for comparison of sensitivities. When the magnitude of one or more columns becomes excessivly large in comparison, the end result is that the resulting model design matrix looks more and more ill-conditioned until the numerics finally fail.
15 Comments
d=dir('median*.mat'); % load the files w/o having to type a zillion characters...
for i=1:numel(d)
load(d(i).name)
end
es=medians_energy_scored_tr; ds=medians_dwi_scored_tr; % shorten names
e=medians_energy_vec.'; d=medians_micro_parameter_vec.';
b=polyfit(d,e,1) % compare polyfit() w/o standardization
polyfit() discovers the problem, too.
[b,~,mu]=polyfit(d,e,1) % and with internally scaled
fitlm(zscore(d),e)
fitlm and polyfit agree with each other with OLS and using standardization on the independent variable as would expect.
histogram(d,20)
While the x values are all distinct, we can observer they are clustered in a couple of regions. It might be better numerically (although a different problem, of course) if they were distributed uniformly...let's see what might happen then...
x=linspace(min(d),max(d),numel(d));
polyfit(x,e,1)
fitlm(x,e)
Well, that's not enough on its own; let's go back to just making the magnitude smaller that we did before...
polyfit(x/1E10,e,1)
and Voila! -- at least the coefficients can again be estimted once the values don't completely swamp precision.
Just one more demonstration...
d=dir('median*.mat'); % load the files w/o having to type a zillion characters...
for i=1:numel(d)
load(d(i).name)
end
%es=medians_energy_scored_tr; ds=medians_dwi_scored_tr; % shorten names
e=medians_energy_vec.'; d=medians_micro_parameter_vec.';
o=ones(size(e));
for n=0:14
x=[d/10^n o]; % the design columns; x term scaled by powers of 10
X=x.'*x;
fprintf('Divisor: %2d, Condition Number: %.4e\n',n, cond(X))
end
x=[zscore(d) o]; % now do the z score to compare
X=x.'*x;
fprintf('Standardized, Condition Number: %.2e\n', cond(X))
Illustrates what happens to the design matrix condition number as the magnitude difference between the constant term column of ones and the x values is reduced. Initially, since we're dividing by a factor of 10, the condition number is reduced by almost exactly a factor of 100 owing to the multiplication; once the difference is down to about 8000 differences in the values begin to appear at the higher level of significant digits besides the pure scaling.
When do use actual z-scoring at the end, the actual use of the data itself to remove the mean and standardize does have the remarkable property of bring the condition number down to almost identically 1.
Manuela
on 1 Jan 2026
x.'*x arises in the matrix form of the O(rdinary)L(east)S(quares) calculation -- hence, in one form or another the solution for the coefficients is dependent upon inverting (or other solution technique) the design matrix. For a Nth order polynomial, the base x is a column-wise matrix
[x.^n x.^(n-1) ... x.^2 x 1]
where x is the vector of the independent variable values and the column of ones is there for the intercept term. (You can solve for a non-intercept model by omitting it). And, to compute the sums of squares, the form is to multiply the transpose by x. In your case of simply a first-order polynomial, you have only [x 1] columns as were used in the examples.
The end result is
X'X B = X'y known as the "normal" equations.
where B is the array of coefficients and y the dependent observations vector. So, to solve for B
inverse(X'X)(X'X) B = inverse(X'X) X'y
I B = inverse(X'X) X'y where I is the identity matrix so
B = inverse(X'X) X'y is the solution for the coefficients, B
That's why the X'X matrix is key and if it is ill-conditioned one has trouble and showed the resulting condition number as a function of the magnitude of x to illustrate.
I don't see it written expressly in the MATLAB doc, but <here's a brief development>; note specifically they arrive at eq (10) at the top of page 3.
Just out of curiosity, why the magnitude of the variable in your example? What does it represent? If it is really of that magnitude, your choices for estimation are likely to be to either transform the units to a set that are much smaller in magnitude or use the standardization technique.
Manuela
on 2 Jan 2026
dpb
on 2 Jan 2026
Not exactly, but similar idea -- if the X values in two columns are identical or exact linear combintions of each other, then the system is identically singular -- the values themselves are all finite but the combination of both makes the system inderminate. In this case, the values aren't actually perfectly correlated, but the large magnitude makes the ability to distinguish between many go away because of loss of precision so it appears more and more like they actually are linearly related as the magnitude gets larger.
As for reporting, you can compute the results using standardized variables but then transform the coefficients back to the original units...it's just algebra.
z = (x-u)/s --> z s = x - u --> x = z s + u
Yes, it's the same model, just using the base variable instead of the standardized.
You could count the cells in terms of 10E9 or 10E12 cells/cu-m instead and see the same t-values.
d=dir('median*.mat'); % load the files w/o having to type a zillion characters...
for i=1:numel(d)
load(d(i).name)
end
%es=medians_energy_scored_tr; ds=medians_dwi_scored_tr; % shorten names
e=medians_energy_vec.'; d=medians_micro_parameter_vec.';
fitlm(d,e)
fitlm(d/1E9,e)
fitlm(d/1E12,e)
fitlm(zscore(d),e)
Do you have a real dataset could attach? One presumes from the plots these are just random values over the expected ranges?
Manuela
on 2 Jan 2026
There's essentially zero correlation of the two variables is why I made that presumption as well as the plots look basically like a pretty-much random scattering of points in the 2D space.
What is the null hypothesis here? An adjusted R-square correlation of only 0.001 is pretty much showing that the two are virtually uncorrelated to each other. The slope is only significant at about the 75% level; it is -0.12+/-0.11 -- the standard error of estimate is virtually the same as the estimate.
Manuela
on 2 Jan 2026
No, but one normally uses regression to predict the response variable from the independent one -- but if the total Rsq shows that the model can explain only about 1% of the total variance, it's certainly going to be of little use for that purpose.
It just doesn't seem like an appropriate analysis given the data unless the null hypothesis is to show they aren't related closely and the slope being virtually statistically indistinguishable from zero makes one wonder just what is the underlying objective here.
IOW, I'm wondering if there might not be a more appropriate statistical approach than regression here, depending upon what the null hypothesis really is. Only you/your lead researcher can answer that, of course, and that is getting far off the forum topic of MATLAB-specific functionality (not tht we aren't quite a ways departed already <vbg> ).
Manuela
on 2 Jan 2026
Manuela
on 2 Jan 2026
In that case I'd suggest corr is the more appropriate. Use the two-output form to get the estimated p-value. This will avoid the issue of the regression calculation matrix inversion of the normal equations entirely.
d=dir('median*.mat'); % load the files w/o having to type a zillion characters...
for i=1:numel(d)
load(d(i).name)
end
%es=medians_energy_scored_tr; ds=medians_dwi_scored_tr; % shorten names
e=medians_energy_vec.'; d=medians_micro_parameter_vec.';
[c,p]=corr([d,e],'Type','Kendall');
[c(2,1) p(2,1)] % corr returns the full correlation matrix, not just the 2-column comparison
I used Kendall's tau instead of default as the previous histogram indicated not normally distributed owing to what appeared to be multiple peaks in the distribution. You'll undoubtedly want to do some distribution testing first to see just what sort of assumptions can be justified.
More Answers (0)
Categories
Find more on Descriptive Statistics in Help Center and File Exchange
