How to find the row wise p values corresponding to each coefficient of linear regression?

I have 3 matrices Y, x1, and x2 each having 3420 rows and 29 columns. I want the row wise linear regression p values of Y on x1, x2, and the interaction between x1 and x2 such that the p value term for each of them is of 3420 rows and 1 column.

I tried the following code:

nsample = 29;
nvar = 3420;
% sample data 
y = rand(nsample, nvar);
x1 = rand(nsample, nvar);
x2 = rand(nsample, nvar);
xInt = x1.*x2;
intercept = ones(nsample, 1);
beta = zeros(nvar, 4);
for i = 1:nvar
    desMat = [intercept, x1(:, i), x2(:, i), xInt(:, i)];
    beta(i, :) = regress(y(:, i), desMat);
end

But beta only gives the linear regression coefficients. I want the p values.

How to get the p values corresponding to the linear regression coefficients of Y on x1, x2, and interaction between x1, x2 as a 3420 rows and 1 column array?

[Jeff Miller]

Do you want the p value of each predictor based on its Type I, II, or III sum of squares (e.g., intro to sums of squares types)?

[Abhishek Chakraborty]

I want Type I sum of squares so that the p value for each predictor is an array containing 3420 rows and 1 column.

Answer 1

[Jeff Miller]

For Type I (i.e., sequential), you have to fit a series of three models, and you need to compute the improvement of each model over the previous one. For this application, I think it is more convenient to use the fitlm command which provides its own intercept & easier-to-use output. So, I think this will do what you are asking for:

nsample = 29;
nvar = 3420;
% sample data 
y = rand(nsample, nvar);
x1 = rand(nsample, nvar);
x2 = rand(nsample, nvar);
xInt = x1.*x2;
% intercept = ones(nsample, 1);  % not needed
beta = zeros(nvar, 4);
type_1_p_value_for_x1 = zeros(nvar, 1);
type_1_p_value_for_x2 = zeros(nvar, 1);
type_1_p_value_for_int = zeros(nvar, 1);

for ivar = 1:nvar
    desMatx1 = [x1(:, ivar)];
    mdlx1 = fitlm(desMatx1,y(:,ivar));
    Fx1 = mdlx1.SSR / mdlx1.MSE;
    type_1_p_value_for_x1(ivar) = 1 - fcdf(Fx1,1,mdlx1.DFE);
    desMatx1x2 = [x1(:, ivar), x2(:, ivar)];
    mdlx1x2 = fitlm(desMatx1x2,y(:,ivar));
    Fx2givenx1 = (mdlx1x2.SSR - mdlx1.SSR) / mdlx1x2.MSE;
    type_1_p_value_for_x2(ivar) = 1 - fcdf(Fx2givenx1,1,mdlx1x2.DFE);
    desMatx1x2int = [x1(:, ivar), x2(:, ivar), xInt(:,ivar)];
    mdlx1x2int = fitlm(desMatx1x2int,y(:,ivar));
    Fintgivenx1x2 = (mdlx1x2int.SSR - mdlx1x2.SSR) / mdlx1x2int.MSE;
    type_1_p_value_for_int(ivar) = 1 - fcdf(Fintgivenx1x2,1,mdlx1x2.DFE);
    beta(ivar, :) = mdlx1x2int.Coefficients.Estimate;
end

[Abhishek Chakraborty]

Thanks a lot, Jeff!