# Which standard deviation to use for p-value of regression coefficients?

Matlab generates regression coefficients for vector time series models ("vgxvarx").

For a given regression coefficient, the p-value is just the area under the t-distribution beyond |t| > (t-statistic), which is easy to do, as per this answer on Matlab Central. The expression is

pValue = 2 * tcdf(-abs(t), dfe).


I have two questions about the conversion from t to p.

The first has to do with the fact that the standard errors are described as maximum likelihood and not corrected to reflect the degrees of freedom. From Wikipedia, it seems that I would need to modify t by a factor sqrt(df/n), where n is the number of observations and df is n less the number of estimated parameters. Can anyone chime in with how advisable this is? My gut says that it boils down to whether I want an estimate of dispersion that is maximum-likelihood or reduced-bias. So I fear that there is probably no right answer, and I'm hoping that someone can dispel that fear.

The second question has to do with the factor of 2 in the above expression for pValue. I think I want the factor of 2 because I'm interested in the area under the t-distribution for both regimes t<-|t-statistic| and t>|t-statistic|. The author of the post weaves in an explanation of 2 standard errors, but I suspect that this is unrelated to seeking the p-value for a specific t-statistic, especially since the value of |t| for 5% significance can depart nontrivially from 2.0 for small sample size. Is this correct?

The above has also been posted to usenet groups matlab and stat.

The data to explore this problem is linked to from this tutorial. The data file is located eggs.txt.

Here is the code that I used to explore the model:

clear variables

fnameDiary = [mfilename '.out.txt'];
if logical(exist(fnameDiary,'file'))
diary off
delete(fnameDiary)
end % if
diary(fnameDiary) % Also turns on diary

if0mean=true % Subtracts out the mean
ifminBias=false % std err normalizes by df rather than n,
% where n=number of observations and df is n less the
% number of parameters estimated (from vgxdisp)
nMaxLag=3 %#ok<NOPTS>
clf

tbChicEgg = readtable( 'eggs.txt', 'Delimiter', '\t' );

seriesNames={'chic','egg'};
varChicEgg = vgxset( 'Series', seriesNames, 'n',2 );
estVarChicEgg(nMaxLag).var = vgxset;

chicEgg = table2array(tbChicEgg(:,seriesNames));
if if0mean
chicEgg = bsxfun( @minus, chicEgg, mean(chicEgg) );
end
chicEgg0 = chicEgg(1:nMaxLag,:);
chicEgg = chicEgg(1+nMaxLag:end,:);
yrs = table2array(tbChicEgg(1+nMaxLag:end,'year'));

subplot(1+nMaxLag,1,1);
plotyy( yrs,chicEgg(:,1) , yrs,chicEgg(:,2) );

for maxLag = 1:nMaxLag

[est.spec est.stdErr est.LLF est.W] = vgxvarx( ...
vgxset( varChicEgg, 'nAR',maxLag ), ...
chicEgg, NaN, chicEgg0, ...
'StdErrType', 'all' ...
);

fprintf('-------------------------\nmaxLag=%g\n',maxLag);
estVarChicEgg(maxLag).var = est %#ok<NOPTS>

subplot(1+nMaxLag,1,1+maxLag)
plotyy(yrs,est.W(:,1),yrs,est.W(:,2))

end

%--------------------------------------------------
diary off
% error('Stopping before return.');
return

• I understand your reasoning, but please do not post the same question to multiple SE sites. – A. Donda May 29 '15 at 23:57
• This is the 2nd time that I've been told this. The 1st time, I looked quickly for a policy and did not find one. I replied that the issue had no consensus, at least in general. This time, I looked harder, and even though opinions are divided, there is a leaning against cross-posting. Furthermore, it is discouraged in the Help. So I will adhere to that. Someone already removed the posting on StackExchange. However, why were the links to usenet removed? – StatSmartWannaB May 30 '15 at 19:38
• Because they add nothing to the question. Another question: did my answer help you? – A. Donda May 30 '15 at 23:39
• But that is the repeated recommendation of all the discussions that I've perused. Also, there is another response on usenet which people may want to see. As for your answer, I had to take the time to investigate alternative explanations. Your part 2 helped....part 1, I've found the answer to, which I will explain in the comments to the answer. I'm not sure how one normally proceeds in such a case. It would be good if people could see both as answers. – StatSmartWannaB May 30 '15 at 23:47
• Maybe you're right. Feel free to re-introduce the links then. – A. Donda May 30 '15 at 23:49

First of all, that answer on Matlab Central is really bad. It doesn't really answer the question, and it gives some bad advice. The formula they give computes the $p$-value for a two-sided test given a $t$-statistic and its degrees of freedom. The information they failed to give is how to obtain the error degrees of freedom dfe from the output of vgxvarx. And saying, don't bother about $p$-values just look at whether $|t| > 2$, ignores a host of problems, among other things how to deal with multiple tests. I'm afraid I cannot clear up all of it, but here's a partial answer.

1) You are right that in estimating the model, error degrees of freedom are lost, and in order to properly compute the $t$-statistic and to properly interpret the $t$-statistic, we would need to know how many error degrees of freedom that are. It appears that using the syntax

vgxdisp(EstSpec, EstStdErrors, 'DoFAdj', true)


you can obtain properly computed $t$-statistics, but the output still does not specify the degrees of freedom, so you still cannot properly interpret the result, and in particular, you cannot compute the $p$-value. Without wading through the internals of vgxvarx and vgxdisp, I don't see a way to determine these degrees of freedom. My recommendation: Ask the MathWorks, and don't accept a crappy answer.

If you know you have a very large number of data points, you can approximately assume that dfe $\rightarrow \infty$, in which case the $t$-distribution becomes the standard normal distribution, and you can approximately compute pValue = 2 * normcdf(-abs(t)).

Regarding the factor you mention, sqrt(df/n), no, this is not for correcting dfs, but to convert from the standard deviation associated with a mean to the standard error of that mean (and regression coefficient estimates can be seen as weighted means). This formula is not appropriate here.

2) Yes, the factor 2 comes from the fact that the $p$-value for a two-sided test is being computed. You are right that this is not related to the "rule of thumb" mentioned there, and you are right that this rule of thumb can be severely off. The rule of thumb is however related to the normal distribution mentioned above: 2 * normcdf(-2) is about 0.0455.

• Thanks, A. Donda. I found evidence that my conjecture in part 1 is (almost!) correct. If you agree with my answer (being composed now and soon posted), could you please modify your part 1? I wish to avoid confusion. If not, we can discuss and clear up what might be behind the discrepancy. P.S. By the way, I have so many questions on the vgx library functions that there is no point trying to contact TMW. – StatSmartWannaB May 31 '15 at 0:16
• I marked this as the acccepted answer, but as a caveat for part #1, please see my answer as well. – StatSmartWannaB May 31 '15 at 23:33

Thanks to A. Donda for the answer to part 2. I was able to investigate my conjecture for part 1 in my original post for this question. Yes, maximum likelihood estimate of the standard error for the coefficient uses a normalization of $n$ (the number of observations) whereas minimum bias estimate uses a normalization of $n$ less the number of coefficients. I counted the number of coefficients from the vgxdisp output for the case of (maximum lag)=3. I then compared the standard errors of the estimated coefficients with the vgxvarx's 'DoFAdj' set to true and to false (thanks, A. Donda), which yields the minimum bias and maximum likelihood estimates, respectively.

My comparison of the two sets of standard errors consisted of multiplying the maximum likelihood standard errors by $n$ to cancel out the normalization by $n$, then dividing by $n$ less the number of parameters (the normalization for minimum bias). Note that I had it backward in my original post.

There is a huge caveat though. I had to increase the number of parameters by exactly 1 in order for the numbers to match. So there must be an extra estimated parameter which isn't listed in the vgxdisp output. Perhaps it is the cross-covariances between the noise sources for the two time series in the model. It only shows up once in the vgxdisp output, but the cross-covariances were considered to be the same parameter because they should be identical. Maybe for regression purposes, they are not.

Here is the code that converts maximum likelihood estimates of standard error to minimum bias estimates, with the number of estimated parameters $nEstParam$ being 1 more than the number listed in vgxdisp's output:

n=74; nEstParam=3*4+3+1; x=sqrt(n/(n-nEstParam)) * ...
[
0.138092
18.6127
0.00102525
0.138188
0.151676
29.8176
0.0011261
0.221377
0.131599
19.4701
0.000977039
0.144553
5.81405e+07
349742
3204.78
]; for i=1:length(x); disp(x(i)), end


P.S. I realize that the solver is reporting poor conditioning of the matrix. I am guessing that it is due to the apparent unit root in the 2nd time series, as hinted by the 1st subplot.

NOTE: As much as this answer confirms what is done in estimating the parameters, it still doesn't answer the question of whether the maximum likelihood or minimum variance estimate is better. Not just better in a general sense (which probably hs no answer), but which is better for estimating p-value. I am not sure whether there is a basis that can be put forth to favour one over the other. Thoughts?

New info: I feel it is important to share that the mystery extra degree of freedom is not due to counting all 4 coefficients in the 2x2 covariance matrix, as suggested by the 3+1 contribution to nEstParam above. This was tested by selecting a diagonal covariance matrix (CovarType='diagonal' for vgxvar). We would expect that replacing 3+1 by 2 in nEstParam would yield the right conversion from ML stderr to min bias. It doesn't, and it turns out that the right conversion factor results from replacing 3+1 by 2+1 in nEstParam above. I suspect that vgxdisp is using the wrong degrees of freedom to calculate stderr because the resulting nEstParam doesn't even match the number of active parameters from vgxcount. I could be wrong, but until the extra DoF can be explained, it's the only possible conclusion. Hence, if I had to extract stderr from vgxdisp, I would request a ML stderr and correct it using the DoF from vgxcount. This is not ideal because it is still guessing, but it's the most sensible guess so far.

• Hi StatSmartWannaB, good that you figured that out. I agree, the t-statistic should be computed using the square root of the unbiased estimator of the variance. I didn't write this explicitly in my answer, but I don't see any disagreement either; 'DofAdj' is the way to go. – A. Donda May 31 '15 at 23:43
• To my understanding, that same number should be the number of degrees of freedom of the resulting test statistics, so I think you should use that number as dfe when computing the p-value. However, I'm not sufficiently familiar with vector autoregressive models to positively say that is the case, and my quick look into the code didn't provide any more insight. What you can do is try to have a closer look yourself and figure out where exactly the listed t-value is being computed, and how. If this agrees with your numerical results, I guess you're good. Sorry I couldn't be of more help. – A. Donda May 31 '15 at 23:47
• The t-value is just the signal to noise (estimated coefficient divided by standard error). I was trying to divine a theoretical or conceptual argument for favouring the standard error from max likelihood or min bias (there's still some bias because taking the square root of the unbiased variance estimate is nonlinear). I think I'll go with minimum bias, because I cannot see why the maximum likelihood is so important, even in a simple generic situation. For example, imagine a probability distribution that is lopsided. Why would you take the peak as the representative value? Anyway, thanks! – StatSmartWannaB Jun 1 '15 at 0:07
• OK, I didn't understand that that is the question. Since it is, there is a clear answer: The standard deviation that enters the $t$-statistic is the square root of the unbiased variance estimate. There is no arguing, that is just the definition. – A. Donda Jun 1 '15 at 0:13
• This definition can be motivated though: Using the ML estimate, the resulting statistic would not have a $t$-distribution, but a scaled $t$-distribution. Or more generally speaking: It would have a distribution that does not solely depend on a single parameter, namely the error degrees of freedom, but also on the total number of samples. – A. Donda Jun 1 '15 at 0:13