comparing regression slopes
Mike Pesavento
I am looking for something similar to a two-sample ttest, of the form:
t= (b1-b2) / (standard error of difference between coefficients)
and all of the various calculations that come with solving for the denominator. This would give me the t value, and I could look up on a table what the corresponding P value is for the given DF. This part is pretty straightforward, but I would rather be able to get a precise p-value for the associated tvalue, and I havent seen a pval calculator in the stats toolbox.
I have also tried using aoctool(), which gives me the slopes and covariance across groups, and then push the stats struct into multcompare(). This last function looks like it is comparing confidence interval overlap to see if slopes are different, and doesnt give a significance associated with it, other than whatever is used to define the confidence intervals.
Long story short, I'm looking for something that can either compare two regression coefficients giving a pval, or find a tool that given a tvalue and df will give me the associated pval. Thanks!
Mike Pesavento
Tom Lane
>statistically different. This seems like a fairly straightforward
>operation, and is present in pretty much ever statistics textbook I've
>found. But I can't seem to find how to do so with the statistics toolbox.
Mike, here are five ways to do it:
1. aoctool coefficients table
load fisheriris
x = meas(51:150,3);
y = meas(51:150,4);
g = species(51:150);
aoctool(x,y,g)
% look at coefficients table, the two group slope coefficients are
% constrained to sum to 0, so just test that one of them is 0 to
% get t=2.08, p=0.0398
2. aoctool anova table
% same as above, but notice the interaction term, which specifies
% a separate slope for each group, has F=4.34, p=0.0398
3. regression using regstats with an X matrix where the coefficient of the
last column represents the difference in slopes between groups 1 and 2, so
just do a t-test on that
group2 = (g=='virginica');
X = [group2, x, x.*group2];
s = regstats(y,X);
[s.tstat.beta, s.tstat.se, s.tstat.t, s.tstat.pval]
% note last row has t=-2.08, p=0.0398
4. regression using regstats, but setting up X in a way that might seem more
natural, with separate slopes for each group, then do an F test for a
difference between them
group1 = ~group2;
X = [group2, x.*group1, x.*group2];
s = regstats(y,X);
[p,F] = linhyptest(s.beta, s.covb, 0, [0 0 1 -1], s.tstat.dfe)
% result is F=4.34, p=0.0398
5. the by-hand way you described
deltab = s.beta(3)-s.beta(4)
se = sqrt(s.covb(3,3) + s.covb(4,4) -2*s.covb(3,4))
t = deltab/se
p = 2*tcdf(-abs(t), s.tstat.dfe)
% result is t=2.08, p=0.0398
Note that for a 1 d.f. test like this, F=t^2.
-- Tom
Mike Pesavento
Tom, thank you for such a thorough clarification of these tests. I think the greatest difficulty I had was in the correct interpretation (and generation!) of the interaction term across all of the tests. Definitely more of a theoretical error on my part. The stats struct that I was getting from aoctool (didnt use the graphical version until your example) as well as from regstats did not have a clear comparison between the slopes.
Part of my confusion comes from the fact that I believe I should be doing multiple regression on my data, with two independent variables (and multiple groups) controlling the result of the dependent variable. In this case, the interaction term would mean something far different, since a given pair of slopes would describe a plane that I would then like to compare to a separate group, obtaining separate estimates of the slopes. From the examples that you gave me, I would assume that it is possible to do higher order interactions between the IVs (in my data, timing and magnitude of input) and the groups to be able to pull out a comparison of the slope along timing for one group and compare it to the slope along timing for the other group.
I have been using regress() to pull out these slopes, and can see how to use regstats to do something similar for a single plane (two IVs, one DV and compare the two slopes of the plane), but not across pairs of slopes between two planes. Is this possible, again by setting up interaction groups between the IVs?
eg:
X=[ group1 x1 x2 x1.*group1 x2.*group1];
s = regstats(y,X);
(similar to example 3 using regstats)
Thank you! You've been a great help in figuring this out so far!
-Mike
Tom Lane
> use regstats to do something similar for a single plane (two IVs, one DV
> and compare the two slopes of the plane), but not across pairs of slopes
> between two planes. Is this possible, again by setting up interaction
> groups between the IVs?
> eg:
> X=[ group1 x1 x2 x1.*group1 x2.*group1];
> s = regstats(y,X);
> (similar to example 3 using regstats)
Mike, I can't say that I completely understand this. But suppose you set up
X this way:
X = [group1 group1.*x1 group2.*x1 group1.*x2 group2.*x2]
Then if you mimic my linhyptest example, you could set up contrast matrices
for a variety of tests:
C1 = [0 0 1 -1 0 0] % each group has the same slope for x1
C2 = [0 0 0 0 1 -1] % each group has the same slope for x2
C3 = [0 0 1 -1 0 0;
0 0 0 0 1 -1] % same x1 slope and same x2 slope across groups
C4 = [0 0 1 0 -1 0;
0 0 0 1 0 -1] % x1 and x2 have the same slope within each group
C5 = [0 0 1 -1 0 0;
0 0 1 0 -1 0;
0 0 1 0 0 -1] % slopes match across predictors and groups
-- Tom
Mike Pesavento
> X this way:
>
> X = [group1 group1.*x1 group2.*x1 group1.*x2 group2.*x2]
>
> Then if you mimic my linhyptest example, you could set up contrast matrices
> for a variety of tests:
YES! Thank you, that was what I needed to pull things together. The masking contrasts within linhyptest() make perfect sense for what I want to do, and I'm sure I will use it many times over in the future. Thank you so much!
-Mike
Akshay
How would you formulate case-4 with covariates?
Best,
Akshay