How to get the standard error of regression coefficient?

[Ivy]

For a linear model Y=XB+u,how can I know the standard error of B if
using
b = regress(Y,X) to estimate B?

Similarly, how can I know the standard error of each element of Bs if
I use
Bs = bootstrp(100,'regress',Y,X)?

Thanks.

[Steven Lord]

Ivy wrote:
> For a linear model Y=XB+u,how can I know the standard error of B if
> using
> b = regress(Y,X) to estimate B?

Do you want confidence intervals for the coefficients B? If so, just
specify more output arguments for your call to REGRESS. The help says, in
part:

REGRESS Multiple linear regression using least squares.
b = REGRESS(y,X) returns the vector of regression coefficients, b,
in the linear model y = Xb, (X is an nxp matrix, y is the nx1
vector of observations).

[B,BINT,R,RINT,STATS] = REGRESS(y,X,alpha) uses the input, ALPHA
to calculate 100(1 - ALPHA) confidence intervals for B and the
residual vector, R, in BINT and RINT respectively. The vector
STATS contains the R-square statistic along with the F and p
values for the regression.

> Similarly, how can I know the standard error of each element of Bs if
> I use
> Bs = bootstrp(100,'regress',Y,X)?

That's a bit more difficult ... you'd probably need to make a backup copy of
bootstrp.m, modify it so that it FEVALs the REGRESS function with multiple
outputs, and formats the extra outputs in the format you need.

--
Steve Lord
sl...@mathworks.com

[Ivy]

I need to know the standard error,not just the confidence interval.
Is there any way to figure it out? Can you show me how to make a

backup copy of bootstrp.m, modify it so that it FEVALs the REGRESS

function with multiple outputs?

In fact my purpose is to construct bootstrap confidence interval.Do
you know any easy way?

Thank you very much.

Ivy

Steven Lord wrote:

>
>
> Ivy wrote:
>> For a linear model Y=XB+u,how can I know the standard error of
B if
>> using
>> b = regress(Y,X) to estimate B?
>

> Do you want confidence intervals for the coefficients B? If so,
just
> specify more output arguments for your call to REGRESS. The help
says, in
> part:
>
> REGRESS Multiple linear regression using least squares.
> b = REGRESS(y,X) returns the vector of regression coefficients, b,
> in the linear model y = Xb, (X is an nxp matrix, y is the nx1
> vector of observations).
>
> [B,BINT,R,RINT,STATS] = REGRESS(y,X,alpha) uses the input, ALPHA
> to calculate 100(1 - ALPHA) confidence intervals for B and the
> residual vector, R, in BINT and RINT respectively. The vector
> STATS contains the R-square statistic along with the F and p
> values for the regression

>> Similarly, how can I know the standard error of each element of
Bs if
>> I use
>> Bs = bootstrp(100,'regress',Y,X)?
>

[Peter Perkins]

> For a linear model Y=XB+u,how can I know the standard error of B if
> using b = regress(Y,X) to estimate B?

REGRESS is set up to return CIs rather than SEs; you probably want to use
soemthing like

stats = regstats(Y,X,'linear',{'beta','covb'}

instead, then take the sqrt of the diagonal of stats.covb. One thing to
remember: 'X' in the call to REGSTATS is not the same as for REGRESS -- it is
the predictor variable only (i.e., a vector), and should not be the design
matrix that includes a column of ones for the intercept term. REGSTATS creates
the design matrix for you.

> Similarly, how can I know the standard error of each element of Bs if
> I use Bs = bootstrp(100,'regress',Y,X)?

Not sure what you mean by "each element of Bs".

If you are trying to get a bootstrap estimate of the SE of the estimated
coefficients in the regression model Y=XB+u (i.e., the same quantity from your
first question), then you really mean "each _column_ of Bs", and you should use
std(Bs).

If you are bootstrapping to try and get at the sampling distribution of those
SEs themselves, then you will need to write a wrapper around REGSTATS (or just
modify REGSTATS itself) to return sqrt(diag(stats.covb)), and run BOOTSTRP on that.

I suspect you want the former, not the latter. Hope this helps.

- Peter Perkins
The MathWorks, Inc.

[Ivy]

I have tried

stats = regstats(Y,X,'linear',{'beta','covb'}

It turned out to be an error message "Too many input arguments." I
checked the syntax in the help it's
"regstats(responses,DATA,'model')".
So I tried
stats = regstats(Y,X,'linear'). It came out another error message "
Too many output arguments." I did use vector for Y and X,not
including a column of ones. I don't know what's wrong.

I do want the distribution of
those SEs themselves. But I don't know how I can modify REGSTATS
itself to return sqrt(diag(stats.covb). Is there any way that I can
call the REGSTATS.m which should be a built-in m-file in matlab so
that I can modify it?

Thank you.

Peter Perkins wrote:
>
>
>> For a linear model Y=XB+u,how can I know the standard error of
B if
>> using b = regress(Y,X) to estimate B?
>

> REGRESS is set up to return CIs rather than SEs; you probably want
to use
> soemthing like
>
> >

> instead, then take the sqrt of the diagonal of stats.covb. One
thing to
> remember: 'X' in the call to REGSTATS is not the same as for
REGRESS --
> it is
> the predictor variable only (i.e., a vector), and should not be the
design
> matrix that includes a column of ones for the intercept term.
REGSTATS
> creates
> the design matrix for you.
>

>> Similarly, how can I know the standard error of each element of
Bs if
>> I use Bs = bootstrp(100,'regress',Y,X)?
>

[Peter Perkins]

> I have tried
> stats = regstats(Y,X,'linear',{'beta','covb'}
> It turned out to be an error message "Too many input arguments." I

>> x = (1:25)';
>> y = 1 + 2*x + randn(size(x));
>>
>> stats = regstats(y,x,'linear',{'beta','covb'})
stats =
source: 'regstats'
beta: [2x1 double]
covb: [2x2 double]
>> stats.beta
ans =
0.5410
2.0389
>> sqrt(diag(stats.covb))
ans =
0.5051
0.0340

> I do want the distribution of
> those SEs themselves. But I don't know how I can modify REGSTATS
> itself to return sqrt(diag(stats.covb). Is there any way that I can
> call the REGSTATS.m which should be a built-in m-file in matlab so
> that I can modify it?

"edit REGSTATS" will allow you to modify it. But putting a wrapper around it is
even simpler. Create an m file containing this:

function outstat = regwrapper(y,x)
stats = regstats(y,x,'linear',{'beta','covb'});
outstat = [stats.beta(:) sqrt(diag(stats.covb))];

and then from the command line,

>> std(bootstrp(1000,'regwrapper',y,x))
ans =
0.4116 0.0247 0.1115 0.0075

The first two values are BS estimates of the sampling standard deviation of the
estimated coefficients. The last two entries are BS estimates of the sampling
standard deviation of the estimated SEs of the (estimated) coefficents.

There is no reason to do any of this with a bootstrap by the way, unless you are
investigating something like robustness to model misspecification. There are
exact distributional results for both the estimated coeffs and their SEs under
the assumption of i.i.d normal errors.

[Ivy]

The reason that I want to get SEs is that I need them to get the
bootstrap confidence interval.The way I know to construct bootstrap
confidence interval is like this:

1. Compute the original estimate of Beta and standard error (Se)
2. For j=1 to B do :
Generate a sample
Compute the estimate of Beta based on the sample: Beta(j)
compute an estimate of the standard deviation of Beta(j): Se(j).
Compute t(j)=(Beta(j)-Beta)/ Se(j)
3. Compute C_{alpha/2} and C_{1-alpha/2} for the sample formed by all
the t's.
4. Construct confidence interval
LO=Beta - C_{1-alpha/2} *Se
HI=Beta - C_{alpha/2} *Se

It's complex.Do you know any better way to construct the bootstrap
confidence interval?

Thanks.

[Peter Perkins]

> The reason that I want to get SEs is that I need them to get the
> bootstrap confidence interval.

Again, under the standard assumptions for linear regression, there are exact
distributions available for the coefficients and their SE's, and so you can
compute CIs without resorting to bootstrap calculations (e.g., use the second
output of regress). So unless you are trying to investigate robustness with
respect to departures from those assumptions, you do not need to bootstrap.

[Ivy]

The problem that I'm trying solve does not satify the standard
assumption for linear regression.The underlying distribution is not
normal.So I need bootstrapping.Do you know any easier way to get the
bootstrap confidence interval?

Thanks.

[Ivy]

It's really strange.I tried the same example you used.

x = (1:25)';
y = 1 + 2*x + randn(size(x));
stats = regstats(y,x,'linear','beta','covb'})

It still came out the error message "Too many input arguments." Is it
because the version of matlab I'm using?My matlab is 6.1.

Peter Perkins wrote:
>
>
>> I have tried
>> stats = regstats(Y,X,'linear',{'beta','covb'}

[Peter Perkins]

> It still came out the error message "Too many input arguments." Is it
> because the version of matlab I'm using?My matlab is 6.1.

Yes, you would need the R13 version of REGSTATS (MATLAB 6.5 and Stats Toolbox
4.0). However, with R12 you can also use GLMFIT to get SEs for a linear model:

>> x = (1:25)';
>> y = 1 + 2*x + randn(25,1);
>>
>> [beta,dev,stats] = glmfit(x,y,'normal');
>> beta
beta =
1.4472e+000
1.9626e+000
>> stats.se
ans =
4.2561e-001
2.8630e-002