pvalues in RLMResults?
Nate Vack
So, I'm playing with robust regression, and notice there's no pvalues
attribute in the RLMResults I get back. There's a tvalues, though.
Apologies for the naivety of this question, but is there any reason I
can't make myself some p's out of those t's?
Also -- thanks for the great package!
-Nate
Skipper Seabold
> Hi all,
>
> So, I'm playing with robust regression, and notice there's no pvalues
> attribute in the RLMResults I get back. There's a tvalues, though.
>
> Apologies for the naivety of this question, but is there any reason I
> can't make myself some p's out of those t's?
>
So, the short answer is that I've never gone back to Huber's book to
figure out how to treat the standard errors to get the "correct" t
values, and whether they are not actually "z" values when based on the
robust covariance matrix.
Please see this brief discussion as well. I will go have a look ASAP
and push a definitive patch and a more complete answer, unless someone
knows better than I.
https://groups.google.com/group/pystatsmodels/browse_thread/thread/b152653115c5967d/d62cc5e4add2b39e
> Also -- thanks for the great package!
> -Nate
You're very welcome. Thanks for the kind words.
Skipper
josef...@gmail.com
> On Fri, Jun 17, 2011 at 1:51 PM, Nate Vack <njv...@wisc.edu> wrote:
>> Hi all,
>>
>> So, I'm playing with robust regression, and notice there's no pvalues
>> attribute in the RLMResults I get back. There's a tvalues, though.
>>
>> Apologies for the naivety of this question, but is there any reason I
>> can't make myself some p's out of those t's?
>>
>
> So, the short answer is that I've never gone back to Huber's book to
> figure out how to treat the standard errors to get the "correct" t
> values, and whether they are not actually "z" values when based on the
> robust covariance matrix.
>
> Please see this brief discussion as well. I will go have a look ASAP
> and push a definitive patch and a more complete answer, unless someone
> knows better than I.
>
> https://groups.google.com/group/pystatsmodels/browse_thread/thread/b152653115c5967d/d62cc5e4add2b39e
R MASS rlm doesn't report p-values either
> res = rlm(stack.loss ~ ., stackloss, psi = psi.bisquare)
> summary(res)
Call: rlm(formula = stack.loss ~ ., data = stackloss, psi = psi.bisquare)
Residuals:
Min 1Q Median 3Q Max
-10.4356 -1.7065 -0.2392 0.8797 6.9326
Coefficients:
Value Std. Error t value
(Intercept) -42.2853 9.5316 -4.4363
Air.Flow 0.9275 0.1081 8.5841
Water.Temp 0.6507 0.2949 2.2068
Acid.Conc. -0.1123 0.1252 -0.8970
Residual standard error: 2.282 on 17 degrees of freedom
sounds a bit strange, any idea why?
Josef
Nate Vack
computed:
http://www.mathworks.com/help/toolbox/stats/robustfit.html
-Nate
josef...@gmail.com
interesting related discussion for Stata, for sandwich not rlm
normal or t-distribution and if t, which dof
http://www.stata.com/support/faqs/stat/sandwich.html
especially the last few paragraphs
'''
The discussion thus far has not focused on what Stata should or could
do, so let’s examine that.
Does StataCorp have a responsibility to protect users from themselves?
'''
Josef
Skipper Seabold
And now that my memory is being jogged, I think SAS's implementation
of L1 regression/M-estimators reports them as Chi-Square distributed.
It looks like Stata rreg reports Z statistics (they used to report t
statistics and use a different covariance matrix many versions
ago...http://www.stata.com/support/faqs/stat/robust.html), but I
haven't reinstalled Stata (and the manuals) yet. Now I remember why I
left it...but I'll try to resolve it or at least make a transparent
decision sooner rather than later.
Skipper
josef...@gmail.com
chisquare sounds strange, that's positive real line only
> It looks like Stata rreg reports Z statistics (they used to report t
> statistics and use a different covariance matrix many versions
> ago...http://www.stata.com/support/faqs/stat/robust.html), but I
> haven't reinstalled Stata (and the manuals) yet. Now I remember why I
> left it...but I'll try to resolve it or at least make a transparent
> decision sooner rather than later.
my impression from other areas, not rlm,
Z (normal distribution) is a theoretically safe choice because it's
asymptotic. For small samples various adjustments work better (for
example in tests) but often the main justification is by Monte Carlo
and not theoretical.
In large samples it doesn't make much difference, in small samples,
the covariance estimate might not be reliable anyway. (
http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-robust-regression.pdf
page 2)
Josef
>
> Skipper
>
Skipper Seabold
All the pvalues for the t-stats are based on |t|, or more generally
|params/bse|.
>> It looks like Stata rreg reports Z statistics (they used to report t
>> statistics and use a different covariance matrix many versions
>> ago...http://www.stata.com/support/faqs/stat/robust.html), but I
>> haven't reinstalled Stata (and the manuals) yet. Now I remember why I
>> left it...but I'll try to resolve it or at least make a transparent
>> decision sooner rather than later.
>
> my impression from other areas, not rlm,
> Z (normal distribution) is a theoretically safe choice because it's
> asymptotic. For small samples various adjustments work better (for
> example in tests) but often the main justification is by Monte Carlo
> and not theoretical.
Right. This is also the justification given often by Allin Cottrell on
the gretl list.
>
> In large samples it doesn't make much difference, in small samples,
> the covariance estimate might not be reliable anyway. (
> http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-robust-regression.pdf
> page 2)
>
Ok, here it goes...the long answer to 'can I make myself some p's out
of those t's'.
I was under the impression that we provide the choices for the robust
covariance to (attempt to) correct for small sample bias. I need to
go back to the Huber text though. This is the paper that I took the
alternative covariance estimators from.
Huber "The 1972 Wald Memorial Lectures: Robust Regression:
Asymptotics, Conjectures, and Monte Carlo"
"There is little hope to build an exact finite sample theory of robust
regression" but Huber finds that the other covariance estimators that
we provide are really close to the monte carlo results with as few as
4 observations per estimated parameter.
Haha, it's been a long time since I went through this paper, but I
recall that it took me several days. And when Peter Huber is using the
words "tedious," "perfectly horrible," and "not easily comparable" to
describe the robust covariance formulae, I fear there is little hope
for me, a mere mortal. Maybe you can make some sense of it, but
nothing jumped out at me.
There is a discussion in Chapter 5 of Cameron and Trivedi, which says
that since asymptotic theory is used for these non-linear extremum
estimators that Wald tests should use chi-square and normal
distributions rather than F- and t- from linear regression. So maybe
SAS reports chi-square in the same sense that t^2 == F in the case of
a hypothesis test of one-parameter in the linear regression framework.
For tests of linear restrictions: "Practictioners frequently instead
use the F-statistic...in the hope that this might provided a better
finite sample approximation."
For tests of a single coefficient: "Formally sqrt(W) [the "t" stat] is
an asymptotic z-statistic, but we use the notation t as it yields the
usual 't-statistic,' the estimate divided by its standard error. In
finite samples, some statistical packages use the standard normal
distribution whereas others use the t-distribution to compute critical
values, p-values, and confidence intervals. Neither is exactly correct
in finite samples, except in the very special case of linear
regression with errors assumed to be normally distributed, in which
case the t-distribution is exact. Both lead to the same results in
infinitely large samples as the t-distribution then collapses to the
standard normal.
So right now, unless anyone objects, to paraphrase StataCorp, I think
we give users just enough rope to hang themselves. I guess we can go
ahead and kick the chair out and say that the parameter estimates are
distributed (standard) normally and use the standard error from the
robust covariance matrix (so overwrite the deprecated? t() method
provided in LikelihoodModelResults). I will push a patch.
Skipper
josef...@gmail.com
Sounds fun :)
A thought: maybe we should rededicate the deprecated t() method to
give t-values and p-values for different options. For example, we
still miss a simple way to get t- and p-values based on the sandwich
estimators.
Then it would be easier for users to choose the flavor (?) of their rope.
The attributes t-values and p-values remain based on the default
choice of the statistics.
Josef
(A sign that a package is not mature is when every question by a user
triggers a design discussion.:)
>
> Skipper
>
Nate Vack
> (A sign that a package is not mature is when every question by a user
> triggers a design discussion.:)
question from a random user triggers a design discussion.
Thanks for all the info, and yeah -- I guess I'll take the rope ;-)
-n
Skipper Seabold
Yes and this is something that has long bothered me. Another option,
that I previously implemented and I think I prefer, is to make mix-in
classes that have t() and/or z() instead of making a possibly
spaghetti code t() at the top-level. That way it requires thinking
from the developer, but the user just sees what is "correct"/made
available. Especially since I've found myself overwriting t() in the
tsa code more frequently.
Skipper
josef...@gmail.com
If t() has more than just the default, it would have to move into
individual result classes, since many things might be model specific
e.g. which sandwiches are available. mixin or not depends on how we
want to reuse the code.
I like t() even for z because I never learned what z should mean, it's
just Student-t(inf), and it's not necessary to check whether a model
defines a t method or a z method.
Josef
>
> Skipper
>
Skipper Seabold
Same problem with a top-level t() right? unless you have something
else in mind. A bunch of options inappropriate for some results
classes and we end up defining a results-specific t() that only has
certain options and has a super call. It's just a matter of where the
re-used code goes. At the top-level base model or in a separate class.
I think the latter has a better chance to be clearer and
self-contained, but I think we're after the same thing.
> I like t() even for z because I never learned what z should mean, it's
> just Student-t(inf), and it's not necessary to check whether a model
> defines a t method or a z method.
>
Z is the standard notation for standard normal and yeah t(n) converges
in distribution to N(0,1) as n -> inf. I find it misleading to call it
t when the distributional assumptions are different. Ie., the p-values
aren't really based on Student's t (using the residual degrees of
freedom). Also, when the summary() says t, and then people say, well
what DOF are you using, I get different p-values... I understand not
wanting to check for t or z, but I think it's a small price to
(*maybe* have to, in some special situation) pay to be explicit and
transparent.
Skipper
josef...@gmail.com
Ok, some of this needs more thought and overview that I have time for
right now (TSA and kids).
maybe a generic method name might then be easier and a place to store
the method, for example LikelihoodModelResults.conf_int infers the
dist from the name of the model, instead, the model could define the
dist as attribute. Besides normal and t, I also have bootstrap (simple
iid case only) in the base.model.ResultsMixin for
GenericLikelihoodModel, but conf_int is not implemented yet for this.
method : string
Not Implemented Yet
Method to estimate the confidence_interval.
"Default" : uses self.bse which is based on inverse Hessian for MLE
"jhj" :
"jac" :
"boot-bse"
"boot_quant"
"profile"
Once we start to implement more bootstrap like Bart is doing, we need
enough expandability. Although this might affect more conf_int then
t=params/bse_of_some_kind. (I don't know yet about bootstrap p-values,
I only looked at them for statistical tests).
I wouldn't mind keeping the pretty output like summary different from
the internal methods, I think our print/string/summary options will
get more verbose over time.
In either case, I feel like a table with models and methods/attributes
will be necessary to get or keep an overview of what is common and
what's different (now that our models don't just inherit all the same
LikelihoodModelResults methods).
Is it possible just to add pvalues to RLMResults for 0.3.0 and
postpone any other redesign for next time?
Josef
>
> Skipper
>
Skipper Seabold
Likewise, minus the kids.
> maybe a generic method name might then be easier and a place to store
> the method, for example LikelihoodModelResults.conf_int infers the
> dist from the name of the model, instead, the model could define the
> dist as attribute. Besides normal and t, I also have bootstrap (simple
This is the approach I took, but it looks like I deleted the branch already.
> iid case only) in the base.model.ResultsMixin for
> GenericLikelihoodModel, but conf_int is not implemented yet for this.
>
> method : string
> Not Implemented Yet
> Method to estimate the confidence_interval.
> "Default" : uses self.bse which is based on inverse Hessian for MLE
> "jhj" :
> "jac" :
> "boot-bse"
> "boot_quant"
> "profile"
>
> Once we start to implement more bootstrap like Bart is doing, we need
> enough expandability. Although this might affect more conf_int then
> t=params/bse_of_some_kind. (I don't know yet about bootstrap p-values,
> I only looked at them for statistical tests).
>
> I wouldn't mind keeping the pretty output like summary different from
> the internal methods, I think our print/string/summary options will
> get more verbose over time.
>
> In either case, I feel like a table with models and methods/attributes
> will be necessary to get or keep an overview of what is common and
> what's different (now that our models don't just inherit all the same
> LikelihoodModelResults methods).
>
> Is it possible just to add pvalues to RLMResults for 0.3.0 and
> postpone any other redesign for next time?
Sure. I don't have any rushed, big refactors in me right now.
Skipper