using mixins for standard errors
Skipper Seabold
re-organization. In particular, I am starting to push some of the
common code in results up to LikelihoodModelResults that we have been
meaning to do for a while - there's still plenty of code reuse to
exploit. Right now, in conf_int in model, we check explicitly for
certain results classes. The check is fragile and new results classes
have to be explicitly added here, which adds an unneeded complication
for new models I think.
https://github.com/statsmodels/statsmodels/blob/master/scikits/statsmodels/base/model.py#L1218
I tried this once before long ago but am opening it up again for what
we want to do. The code is setup to check an attribute confint_dist,
which we can make sure each model has. This has the shortfall, that
all parameter test statistics are called tvalues whether they are or
not. I made some mixin classes though to be explicit. If theory
dictates that the errors are asymptotically normal, we provide
zvalues. Else, we provide tvalues. This is quite often a subjective
judgement, but I've tended to follow common practices here and the
general rule of thumb MLE estimators -> zvalues, regression and
approximate MLE estimators (e.g., conditional sum of squares) ->
tvalues. The mixin classes are simply
class NormalErrors(object):
@cache_readonly
def zvalues(self):
return self.params / self.bse
class TErrors(object):
@cache_readonly
def tvalues(self):
return self.params / self.bse
The conf_int only needs to check if the Results instance has a zvalues
or tvalues attribute (could be extended to other distributions, though
naming might become an issue, chi2values?). This seems cleaner and the
way to go to me, but what am I missing? It already works well for the
discrete changes.
Aside: I am now thinking that Exact MLE and CSS in ARIMA should be
split up via an intermediate class similar to what I'm doing with
discrete. It's already much cleaner and we can get rid of a lot of the
if 'method' checks.
Skipper
Skipper Seabold
> The mixin classes are simply
>
> class NormalErrors(object):
> @cache_readonly
> def zvalues(self):
> return self.params / self.bse
>
> class TErrors(object):
> @cache_readonly
> def tvalues(self):
> return self.params / self.bse
>
Obviously pvalues can be added here as well or could be pushed up to
LikelihoodModel and check the attribute much like the confidence
intervals.
Skipper
Wes McKinney
I'm fine with this sort of pattern. You could also set class
attributes instead of having a list to check like in the referenced
code-- but that doesn't address the tvalues/zvalues thing which is
cleaner.
Statistical aside, why is self.params / self.bse said to be normally
distributed in the RLM case?
josef...@gmail.com
Two issues
I prefer tvalues to zvalues
I'm in favor of moving things out of LikelihoodModelResults into mixin classes.
to the second:
VAR is missing conf_int and I looked at how MNLogit defines those,
essentially call super (i.e. LikelihoodModelResults) and reshape.
Mixin classes would make it easier to share some common code across
models that do not all subclass LikelihoodModelResults.
LikelihoodModelResults currently does both handle MLE and work as
multivariate t/normal distribution for the estimated parameters.
I think the multivariate t/normal distribution handling should go out
into a mixin class.
normal versus t
during the summary() rewrite I looked at this, and I think I prefer
tvalues now across the board. The summary helper function takes as an
argument whether the values are t or normal , "use_t=True" for the
header naming.
I think I sent an email in favor of tvalues.
Only partially checked: tvalues might in many cases have better small
sample performance then zvalues.
Even in the linear_model case, as soon as we include lagged endogenous
variables or have other deviations from the basic regression model
with normal errors, we are in asymptotic land. Then, from the
theoretical basis, we should only use the normal distribution. But if
nobs is large (comments in code I've seen >100), then there is no
practical difference between normal and t.
normal = t(df=inf)
practical: LikelihoodModelResults assumes df_resid is defined, and
essentially assumes t-distribution, e.g. for f_test and t_test, and
pvalues. Currently only conf_int checks for normal versus t
the definition of tvalues is the same (it's just the naming that is
different), it's the distribution that is different. I think it's
easier to keep an attribute (or even df_resid=np.inf , not in favor)
that can be used in all methods that interpret the distribution of
params. Especially, if specifying the default distribution use_t=True
or use_normal is defined in all models, then we don't need the hasattr
check.
If we really want to separate t and normal, then LikelihoodModel
should get either an optional normal or t mixin, however ...
Another question, if we use a sandwich estimator in conf_int, pvalues,
t-test, f_test, should we also use the normal or keep using the
t-distribution (in OLS or similar).
For the choice of whether to use t or z, I would prefer if it can be
made into an option for the models, rather then hardcoding the choice
by using either one mixin.
>
> Aside: I am now thinking that Exact MLE and CSS in ARIMA should be
> split up via an intermediate class similar to what I'm doing with
> discrete. It's already much cleaner and we can get rid of a lot of the
> if 'method' checks.
two classes or just two different methods? I don't remember the details.
Josef
>
> Skipper
>
Skipper Seabold
> On Thu, Nov 3, 2011 at 12:47 PM, Skipper Seabold <jsse...@gmail.com> wrote:
>> On Thu, Nov 3, 2011 at 12:34 PM, Skipper Seabold <jsse...@gmail.com> wrote:
>>> The mixin classes are simply
>>>
>>> class NormalErrors(object):
>>> @cache_readonly
>>> def zvalues(self):
>>> return self.params / self.bse
>>>
>>> class TErrors(object):
>>> @cache_readonly
>>> def tvalues(self):
>>> return self.params / self.bse
>>>
>>
>> Obviously pvalues can be added here as well or could be pushed up to
>> LikelihoodModel and check the attribute much like the confidence
>> intervals.
>>
>> Skipper
>>
>
> I'm fine with this sort of pattern. You could also set class
> attributes instead of having a list to check like in the referenced
> code-- but that doesn't address the tvalues/zvalues thing which is
> cleaner.
>
Yeah, that's what I'm thinking right now as well.
> Statistical aside, why is self.params / self.bse said to be normally
> distributed in the RLM case?
>
Answered here: http://groups.google.com/group/pystatsmodels/msg/c2953877539e6bb5
Skipper
Skipper Seabold
Sure. Yesterday I made a backwards incompatible change to MNLogit to
list equations in columns rather than rows so that it's like VAR. This
will allow us to re-use code for multiple equations models. Slowly
working in this direction, generic TimeSeriesModel was first on my
TODO list. Caveat, this won't work for multi-equation models that
don't share exog - ie., Heckman two-step, Simulation Equation Models,
etc. We'll need to use a list or a dict. Probably a list, so at least
the indexing can be similar, though I can be convinced otherwise.
> normal versus t
>
> during the summary() rewrite I looked at this, and I think I prefer
> tvalues now across the board. The summary helper function takes as an
> argument whether the values are t or normal , "use_t=True" for the
> header naming.
Sure, but these aren't mutually exclusive, and use_t can be
determined programmatically - might give more code re-use.
> I think I sent an email in favor of tvalues.
> Only partially checked: tvalues might in many cases have better small
> sample performance then zvalues.
Or it might be theoretically the incorrect distribution.
>
> Even in the linear_model case, as soon as we include lagged endogenous
> variables or have other deviations from the basic regression model
> with normal errors, we are in asymptotic land. Then, from the
> theoretical basis, we should only use the normal distribution. But if
> nobs is large (comments in code I've seen >100), then there is no
> practical difference between normal and t.
>
> normal = t(df=inf)
Sure. My understanding is that - don't have the reference now - but I
mention it in the p-values RLM thread, is that linear regression with
normally distributed errors -> t-distribution is the correct
distribution.
>
> practical: LikelihoodModelResults assumes df_resid is defined, and
> essentially assumes t-distribution, e.g. for f_test and t_test, and
> pvalues. Currently only conf_int checks for normal versus t
I think f_test, t_test should become mix-in classes as well, where the
Wald test uses either the F distribution or the Chi-square
distribution depending on theory. We also don't really provide LR or
LM tests either yet I don't think, but these should be mix-ins as
well. We can provide a mix-in aggregator class so we don't have
results classes inheriting from all over the place.
> the definition of tvalues is the same (it's just the naming that is
> different), it's the distribution that is different. I think it's
> easier to keep an attribute (or even df_resid=np.inf , not in favor)
> that can be used in all methods that interpret the distribution of
> params. Especially, if specifying the default distribution use_t=True
> or use_normal is defined in all models, then we don't need the hasattr
> check.
Splitting hairs here, but if you code a new model and you forget to
set the attribute, then you get an answer that may be slightly off
without realizing it (tests!), but if you're explicitly forced to
include an errors mix-in class then this will never come up. We could
even push conf_int to the mix-in classes so there's not hasattr check.
I'm not convinced that use_t is preferable to just being explicit and
naming them tvalues of zvalues.
> If we really want to separate t and normal, then LikelihoodModel
> should get either an optional normal or t mixin, however ...
>
> Another question, if we use a sandwich estimator in conf_int, pvalues,
> t-test, f_test, should we also use the normal or keep using the
> t-distribution (in OLS or similar).
>
I'll have to check, but I don't see why this changes the distribution
of the parameters. IIUC, it just corrects the variance estimate for
possible bias.
> For the choice of whether to use t or z, I would prefer if it can be
> made into an option for the models, rather then hardcoding the choice
> by using either one mixin.
>
Give the user more rope? Why?
>>
>> Aside: I am now thinking that Exact MLE and CSS in ARIMA should be
>> split up via an intermediate class similar to what I'm doing with
>> discrete. It's already much cleaner and we can get rid of a lot of the
>> if 'method' checks.
>
> two classes or just two different methods? I don't remember the details.
>
Two classes with a single parent class. Only need to overwrite a few
methods and have a few additional ones in each. See the discrete
refactor later today for how much nicer it looks, but I have to run
for the moment.
Skipper
josef...@gmail.com
normally distributed errors and exogenous variables, I think only
asymptotic normal results are available for FGLS and for stochastic
regressors (maybe some of the latter cases can still have a
conditional exact normal distribution).
I'm starting to get confused
why should beta distibuted asy normal prevent us from defining tvalues
beta/bse and use the t-distriubtion as approximation to the asymptotic
distribution of the t-value?
If we only want to use the normal then the hypothesis would also be
directly on beta.
>
>>
>> practical: LikelihoodModelResults assumes df_resid is defined, and
>> essentially assumes t-distribution, e.g. for f_test and t_test, and
>> pvalues. Currently only conf_int checks for normal versus t
>
> I think f_test, t_test should become mix-in classes as well, where the
> Wald test uses either the F distribution or the Chi-square
> distribution depending on theory. We also don't really provide LR or
> LM tests either yet I don't think, but these should be mix-ins as
> well. We can provide a mix-in aggregator class so we don't have
> results classes inheriting from all over the place.
compare_LR is currently in linear_model but should go into the
LikelihoodModelResultslevel.
There are some specific LM tests in the specification tests, but I
don't remember a general formula or specifically for linear
restrictions.
>
>> the definition of tvalues is the same (it's just the naming that is
>> different), it's the distribution that is different. I think it's
>> easier to keep an attribute (or even df_resid=np.inf , not in favor)
>> that can be used in all methods that interpret the distribution of
>> params. Especially, if specifying the default distribution use_t=True
>> or use_normal is defined in all models, then we don't need the hasattr
>> check.
>
> Splitting hairs here, but if you code a new model and you forget to
> set the attribute, then you get an answer that may be slightly off
> without realizing it (tests!), but if you're explicitly forced to
> include an errors mix-in class then this will never come up. We could
> even push conf_int to the mix-in classes so there's not hasattr check.
I'm not sure whether slightly off in this case has a clear
interpretation. We might be slightly off compared to *some* other
packages, but it would be just a different approximation for the small
sample case.
My impression is that this story is more a documentation issue, that
we have clearly specified what is used in the different cases to avoid
having to figure out slightly divergent results across packages, than
a question of clear statistical right or wrong results.
>
> I'm not convinced that use_t is preferable to just being explicit and
> naming them tvalues of zvalues.
>
>> If we really want to separate t and normal, then LikelihoodModel
>> should get either an optional normal or t mixin, however ...
>>
>> Another question, if we use a sandwich estimator in conf_int, pvalues,
>> t-test, f_test, should we also use the normal or keep using the
>> t-distribution (in OLS or similar).
>>
>
> I'll have to check, but I don't see why this changes the distribution
> of the parameters. IIUC, it just corrects the variance estimate for
> possible bias.
for example, Greene 5th edition, equation 11-9 and end of section
11.5, equ. 17, it's only an Est.Asy.Var
>
>> For the choice of whether to use t or z, I would prefer if it can be
>> made into an option for the models, rather then hardcoding the choice
>> by using either one mixin.
>>
>
> Give the user more rope? Why?
example: non-linear regression inherits most of linear_model results.
Greene p.176 on testing for linear restriction on the parameters
"...., so the F distribution is only approximate"
"It should be noted that the small-sample behavior of W <JP:Wald
statistic> can be erratic, and the more conservative F statistic might
be preferable if the sample is not to large"
Josef
>
>>>
>>> Aside: I am now thinking that Exact MLE and CSS in ARIMA should be
>>> split up via an intermediate class similar to what I'm doing with
>>> discrete. It's already much cleaner and we can get rid of a lot of the
>>> if 'method' checks.
>>
>> two classes or just two different methods? I don't remember the details.
>>
>
> Two classes with a single parent class. Only need to overwrite a few
> methods and have a few additional ones in each. See the discrete
> refactor later today for how much nicer it looks, but I have to run
> for the moment.
The discrete case refactoring is also good because they are really
different models, CSS and MLE are just two different estimators for
the same model.
Reducing case dependent sequences of `if`s makes clearly a more readable design.
Josef
>
> Skipper
>
Skipper Seabold
I'm not sure I see your point. I think we should distinguish between t
and normal. Purity (when available in this case) over practicality.
>
>
>>
>>>
>>> practical: LikelihoodModelResults assumes df_resid is defined, and
>>> essentially assumes t-distribution, e.g. for f_test and t_test, and
>>> pvalues. Currently only conf_int checks for normal versus t
>>
>> I think f_test, t_test should become mix-in classes as well, where the
>> Wald test uses either the F distribution or the Chi-square
>> distribution depending on theory. We also don't really provide LR or
>> LM tests either yet I don't think, but these should be mix-ins as
>> well. We can provide a mix-in aggregator class so we don't have
>> results classes inheriting from all over the place.
>
> compare_LR is currently in linear_model but should go into the
> LikelihoodModelResultslevel.
>
> There are some specific LM tests in the specification tests, but I
> don't remember a general formula or specifically for linear
> restrictions.
>
Ah, good. Yes, this should be made general.
>>
>>> the definition of tvalues is the same (it's just the naming that is
>>> different), it's the distribution that is different. I think it's
>>> easier to keep an attribute (or even df_resid=np.inf , not in favor)
>>> that can be used in all methods that interpret the distribution of
>>> params. Especially, if specifying the default distribution use_t=True
>>> or use_normal is defined in all models, then we don't need the hasattr
>>> check.
>>
>> Splitting hairs here, but if you code a new model and you forget to
>> set the attribute, then you get an answer that may be slightly off
>> without realizing it (tests!), but if you're explicitly forced to
>> include an errors mix-in class then this will never come up. We could
>> even push conf_int to the mix-in classes so there's not hasattr check.
>
> I'm not sure whether slightly off in this case has a clear
> interpretation. We might be slightly off compared to *some* other
> packages, but it would be just a different approximation for the small
> sample case.
>
I don't really recall too many cases where t vs normal wasn't
consistent across packages. At least with the bread and butter
estimators.
> My impression is that this story is more a documentation issue, that
> we have clearly specified what is used in the different cases to avoid
> having to figure out slightly divergent results across packages, than
> a question of clear statistical right or wrong results.
>
Justification given by Allin Cottrell on why they changed from
t-statistics to z-statistics for ARMA across versions.
"We decided it was best to display these as z-statistics for
estimators whose validity is based on asymptotic theory and where
standard errors are based on ML estimation. That may be debatable,
but it's a legitimate choice."
>>
>> I'm not convinced that use_t is preferable to just being explicit and
>> naming them tvalues of zvalues.
>>
>>> If we really want to separate t and normal, then LikelihoodModel
>>> should get either an optional normal or t mixin, however ...
>>>
>>> Another question, if we use a sandwich estimator in conf_int, pvalues,
>>> t-test, f_test, should we also use the normal or keep using the
>>> t-distribution (in OLS or similar).
>>>
>>
>> I'll have to check, but I don't see why this changes the distribution
>> of the parameters. IIUC, it just corrects the variance estimate for
>> possible bias.
>
> for example, Greene 5th edition, equation 11-9 and end of section
> 11.5, equ. 17, it's only an Est.Asy.Var
>
>>
>>> For the choice of whether to use t or z, I would prefer if it can be
>>> made into an option for the models, rather then hardcoding the choice
>>> by using either one mixin.
>>>
>>
>> Give the user more rope? Why?
>
> example: non-linear regression inherits most of linear_model results.
>
Not sure what you're after here. If they're mix-ins, then you can
choose what it inherits. I've always seen t-ratios and p-values and
chi2 in post-estimation tests for NLLS anyway, but read on.
> Greene p.176 on testing for linear restriction on the parameters
>
> "...., so the F distribution is only approximate"
>
> "It should be noted that the small-sample behavior of W <JP:Wald
> statistic> can be erratic, and the more conservative F statistic might
> be preferable if the sample is not to large"
>
I don't have the reference in front of me, but sure. This is the
somewhat similar issue with RLM. I think this is another who's your
target user. If you're debating about the distribution of your test
statistics, you can call the distribution yourself on the returned
test statistic, but the proliferation of options might just serve to
confuse the average user. Stata for instance, always gives you a Chi2
p-value after lrtest.
To be honest, I'm not even really sure what we're debating at this
point, and I'd really just rather make a decision and get back to
coding. So let's stick to the practical. From where I sit, the
practical debate is
1) mix-ins provide sensible and explicit defaults that the average
user will be comfortable with. The test statistic is there if you want
to use another distribution. And I think it's a 'cleaner' solution in
that it gives you less 'attributes that models have to have' that you
have to remember, which lowers the barrier for other developers who
don't know all the machinery always and everywhere (aside: this is why
I'm not entirely happy with the wrapper implementation - you have to
hand code the directions to the wrapper for each results class, though
I haven't come up with a different solution).
2) OTOH, code reuse that relies on the test statistics would have to
be written in such a way as to be agnostic between tvalue, zvalue,
etc. This would mean, either attaching a distribution flag (use_t),
attach the actual distribution confint_dist, or a get_test_stat helper
function that gets the correct attribute, which isn't really all that
bad.
Anything I'm missing?
The theoretical debate is not going to be solved and that's why I
proposed the (rule of thumb) distinction in the original post. Making
a decision on this front doesn't keep power users from doing something
else.
Skipper
Skipper Seabold
> different models, CSS and MLE are just two different estimators for
> the same model.
> Reducing case dependent sequences of `if`s makes clearly a more readable design.
>
It'll still appear as one estimator to the users, but if it'll work I
think the code separation will be nice.
josef...@gmail.com
What I'm arguing is that if you want purity, then there is almost no
case in economics or most other fields outside of experiments with
fixed regressors where we could use the t-distribution.
t-ratio, tvalues params/bse is mostly a t-distribution concept,
otherwise we could just use params[i] is stats.norm.cdf(beta_,
scale=np.sqrt(bse) (ok we can standardize)
(just being picky and splitting hairs)
emphasis on *a* legitimate choice
I don't like making it conditional on the choice MLE versus Least
Squares. If I do ARMA with conditional leastsquares, I use t, with
conditional MLE I use normal distribution even though they are
asymptotically equivalent?
However, users and I might also prefer better small sample behavior to
asymptotic theory.
>
>>>
>>> I'm not convinced that use_t is preferable to just being explicit and
>>> naming them tvalues of zvalues.
>>>
>>>> If we really want to separate t and normal, then LikelihoodModel
>>>> should get either an optional normal or t mixin, however ...
>>>>
>>>> Another question, if we use a sandwich estimator in conf_int, pvalues,
>>>> t-test, f_test, should we also use the normal or keep using the
>>>> t-distribution (in OLS or similar).
>>>>
>>>
>>> I'll have to check, but I don't see why this changes the distribution
>>> of the parameters. IIUC, it just corrects the variance estimate for
>>> possible bias.
>>
>> for example, Greene 5th edition, equation 11-9 and end of section
>> 11.5, equ. 17, it's only an Est.Asy.Var
>>
>>>
>>>> For the choice of whether to use t or z, I would prefer if it can be
>>>> made into an option for the models, rather then hardcoding the choice
>>>> by using either one mixin.
>>>>
>>>
>>> Give the user more rope? Why?
>>
>> example: non-linear regression inherits most of linear_model results.
>>
>
> Not sure what you're after here. If they're mix-ins, then you can
> choose what it inherits. I've always seen t-ratios and p-values and
> chi2 in post-estimation tests for NLLS anyway, but read on.
which pvalues, t or normal ?
see at end the subclass might inherit attributes without overwriting them.
But I'm not sure whether I want t or normal distribution, I used t
distribution so far.
>
>> Greene p.176 on testing for linear restriction on the parameters
>>
>> "...., so the F distribution is only approximate"
>>
>> "It should be noted that the small-sample behavior of W <JP:Wald
>> statistic> can be erratic, and the more conservative F statistic might
>> be preferable if the sample is not to large"
>>
>
> I don't have the reference in front of me, but sure. This is the
> somewhat similar issue with RLM. I think this is another who's your
> target user. If you're debating about the distribution of your test
> statistics, you can call the distribution yourself on the returned
> test statistic, but the proliferation of options might just serve to
> confuse the average user. Stata for instance, always gives you a Chi2
> p-value after lrtest.
>
> To be honest, I'm not even really sure what we're debating at this
> point, and I'd really just rather make a decision and get back to
> coding. So let's stick to the practical. From where I sit, the
> practical debate is
I think, the main debatable question is whether it's "either-or" (in
what's implemented) or what's the default.
I'm in favor of mix-ins, but what's in the mixin? and if there are
several, are they exclusive?
Instead of providing just "one solution" (Zen), I prefer to go in the
direction of stata and add 500 (?) options (sandwiches, bootstrap,
...), but defining defaults that work for standard data analysis.
I like now summary() as the main tool to provide "interpretation help"
and providing default choices and diagnostic analysis if a user is
less experienced with or less opinionated about the statistical
details.
For example, I wouldn't like it if your mix-ins take away t_test and
f_test and only provide Wald test and similar for some models.
The alternative to mixins and other methods in the results, would be
to outsource more into standalone functions, that could be called with
any result instance or relevant attributes of it. (I saw Stata session
protocols that use post-estimation analysis because the defaults isn't
what the user wanted, and it was only one additional command.)
Some results like confidence intervals are not so easy to quickly
implement even as a power user.
For developers I think it is important to have tested mixins, classes
that can be inherited, or functions that can be called and attached to
a method. Something that I tried to provide with
GenericLikelihoodModel, and something where the intermediate discrete
level will help a lot. (I can dig out my conditional multinomial logit
and inherit from the appropriate classes.)
>
> 1) mix-ins provide sensible and explicit defaults that the average
> user will be comfortable with. The test statistic is there if you want
> to use another distribution. And I think it's a 'cleaner' solution in
> that it gives you less 'attributes that models have to have' that you
> have to remember, which lowers the barrier for other developers who
> don't know all the machinery always and everywhere (aside: this is why
> I'm not entirely happy with the wrapper implementation - you have to
> hand code the directions to the wrapper for each results class, though
> I haven't come up with a different solution).
>
> 2) OTOH, code reuse that relies on the test statistics would have to
> be written in such a way as to be agnostic between tvalue, zvalue,
> etc. This would mean, either attaching a distribution flag (use_t),
> attach the actual distribution confint_dist, or a get_test_stat helper
> function that gets the correct attribute, which isn't really all that
> bad.
I think a flag as attribute will also be useful because it can be
immediately used in any (new) method that is written.
One issue on inheritance
If a result class defines either zvalues or tvalues, then any subclass
that uses the other mixin still inherits the attribute/method from the
superclass(es), so in this case it would have both tvalues and zvalues
(unless we try del attribute in the subclass to remove superclass
attributes).
A hasattr check wouldn't be conclusive in this case.
(This might be the case if a non-linear version with use_t=False
inherits from a linear version with use_t=True.)
details:
I don't like zvalues as a term, it doesn't exist in econometrics (?,
at least I never heard it before statsmodels or scipy.stats), Greene
doesn't have an index entry for it, and I see t-ratio used whether
it's t or normal distribution.
On the other hand, pvalues is a very ambiguous name without qualifier,
I wouldn't mind replacing it with tpvalues and zpvalues or better
pvalues_t, pvalues_norm or "params_pvalues_t" "params_pvalues_norm" to
have it fully qualified :)
(aside: do we ever switch away from bse, or is it too traditional? If
we want to switch to params_se or something like that, then we should
wait too long.)
Sorry for being so longwinded, but I just don't think replacing
tvalues by zvalues for some models solves our problems with this.
Josef
josef...@gmail.com
both GLM and Discrete used the t-distribution in pvalues, and because
I only looked at those, also in summary(),
only in conf_int the normal distribution is used.
RLM has normal distribution throughout.
Josef
Skipper Seabold
Basically, I agree with what you're saying, and I'm willing to agree
that replacing tvalues with zvalues is unnecessary, and we should
retain use_t in summary() for printing results and interpretation.
This probably isn't a case for mix-ins.
Echoing what you've already said (and for my own edification), AFAICT
the t-distribution is never used outside of the linear regression
model. Open to being proven wrong. And it's only appropriate as the
exact distribution for the t-statistics in the linear regression model
iff the errors are normally distributed. A situation Stock and Watson
say is "rarely applicable in economics." Instead, the normal
approximation should be used based on an appeal to the central limit
theorem. Also, if the errors are not normally distributed, then the
OLS estimator is less efficient than MLE and we lose the normally
distributed errors -> t-distribution, and should appeal to CLT, ie.,
using the standard normal. The discussions in Stock and Watson and
Greene say that the t-distribution is used in software packages for
the linear model due to history and out of a notion of conservatism in
small samples (if we aren't absolutely sure about the normality of the
errors), respectively.
This leads me to think that pvalues and conf_int in
LikelihoodModelResults should all use the normal distribution, and we
should overwrite these methods in RegressionResults. Can we agree on
this? I need to go back and replace the use of t in some of the other
models (it's still incorrectly t in ARMA in addition to where you've
pointed out)
I think this is the way to go, similar to your information criteria
additions. Have the functions separate, and pull them into a class for
methods. A lot of methods could be used on their own as functions, and
this would also provide another entry-level for users who just want
them.
I'd argue here that we may not have gotten the inheritance right if
this is the case, for example
Skeleton Class
| |
| Skeleton Class + MixIns = SomeResults
|
Skeleton Class + OtherMixins = OtherResults
> A hasattr check wouldn't be conclusive in this case.
> (This might be the case if a non-linear version with use_t=False
> inherits from a linear version with use_t=True.)
>
> details:
> I don't like zvalues as a term, it doesn't exist in econometrics (?,
> at least I never heard it before statsmodels or scipy.stats), Greene
> doesn't have an index entry for it, and I see t-ratio used whether
> it's t or normal distribution.
>
I guess I heard it in a statistics class first for standard normal
variables. http://en.wikipedia.org/wiki/Standard_score
"Standard scores are also called z-values, z-scores, normal scores,
and standardized variables; the use of "Z" is because the normal
distribution is also known as the "Z distribution"."
Stata mentions t vs. z in the docs here:
http://www.stata.com/support/faqs/stat/oneside.html
Greene talks about z in the section on confidence intervals. 7th
edition, section 4.5.1, Forming A Confidence Interval for a
Coefficient
> On the other hand, pvalues is a very ambiguous name without qualifier,
> I wouldn't mind replacing it with tpvalues and zpvalues or better
> pvalues_t, pvalues_norm or "params_pvalues_t" "params_pvalues_norm" to
> have it fully qualified :)
>
> (aside: do we ever switch away from bse, or is it too traditional? If
> we want to switch to params_se or something like that, then we should
> wait too long.)
>
I think we should get away from bse, though this is going to require
deprecation this late in the game.
> Sorry for being so longwinded, but I just don't think replacing
> tvalues by zvalues for some models solves our problems with this.
>
I agree (with both sentiments).
Skipper
josef...@gmail.com
> Ok, did some homework and am back on the machine with these unpushed changes.
>
> Basically, I agree with what you're saying, and I'm willing to agree
> that replacing tvalues with zvalues is unnecessary, and we should
> retain use_t in summary() for printing results and interpretation.
> This probably isn't a case for mix-ins.
>
> Echoing what you've already said (and for my own edification), AFAICT
> the t-distribution is never used outside of the linear regression
> model. Open to being proven wrong. And it's only appropriate as the
> exact distribution for the t-statistics in the linear regression model
> iff the errors are normally distributed.
and we can condition on the exog, which rules out lagged endogenous,
AR (with OLS/GLS), even if the errors are normal. (IIRC)
> A situation Stock and Watson
> say is "rarely applicable in economics." Instead, the normal
> approximation should be used based on an appeal to the central limit
> theorem. Also, if the errors are not normally distributed, then the
> OLS estimator is less efficient than MLE
if MLE is based on the correctly specified likelihood, OLS has weaker
assumptions (as does quasi-MLE)
and we lose the normally
> distributed errors -> t-distribution, and should appeal to CLT, ie.,
> using the standard normal. The discussions in Stock and Watson and
> Greene say that the t-distribution is used in software packages for
> the linear model due to history and out of a notion of conservatism in
> small samples (if we aren't absolutely sure about the normality of the
> errors), respectively.
I'm not sure I agree with the conclusion, that given that we have only
asymptotic results with normality based on CLT, that we should use the
normal distribution in finite samples, t-distribution is also an
approximation that converges to the normal as df goes to infinity.
I think what to use should be more a question of what's the better
small sample approximation. In scipy.stats I got used to having lot's
of small sample correction factors.
The debate is a bit like what should be the default ddof for numpy.var.
I don't really care much, and I'm fine with using the normal
distribution as default everywhere, but I like the option to change
use_t=True
(which also has the advantage that I could do a quick MonteCarlo with
a given design matrix and see which is better.)
Using the normal by default, but also be a better signal to users that
we don't have an exact small sample distribution.
I never checked whether it would make much of a difference except for
very small samples.
(bootstrap or other approximation might be more important in the small
sample case anyway.)
>
> This leads me to think that pvalues and conf_int in
> LikelihoodModelResults should all use the normal distribution, and we
> should overwrite these methods in RegressionResults. Can we agree on
> this? I need to go back and replace the use of t in some of the other
> models (it's still incorrectly t in ARMA in addition to where you've
> pointed out)
I still think the LikelihodModelResults conf_int should base the
results on the use_t or use_norm if we want to change the name, we can
have a generic top-level implementation, and it can be changed case by
case.
compare to conf_int, pvalues are trivial to calculate, but an if
use_t/use_normal would make the results more consistent.
I don't know yet, in the non-linear case that subclasses the linear
case, I only needed to overwrite a few methods in the result class, so
putting an intermediate level might not so necessary. Result classes
with two methods the rest is inherited? possible but not necessary.
>
>> A hasattr check wouldn't be conclusive in this case.
>> (This might be the case if a non-linear version with use_t=False
>> inherits from a linear version with use_t=True.)
>>
>> details:
>> I don't like zvalues as a term, it doesn't exist in econometrics (?,
>> at least I never heard it before statsmodels or scipy.stats), Greene
>> doesn't have an index entry for it, and I see t-ratio used whether
>> it's t or normal distribution.
>>
>
> I guess I heard it in a statistics class first for standard normal
> variables. http://en.wikipedia.org/wiki/Standard_score
>
> "Standard scores are also called z-values, z-scores, normal scores,
> and standardized variables; the use of "Z" is because the normal
> distribution is also known as the "Z distribution"."
>
> Stata mentions t vs. z in the docs here:
> http://www.stata.com/support/faqs/stat/oneside.html
>
> Greene talks about z in the section on confidence intervals. 7th
> edition, section 4.5.1, Forming A Confidence Interval for a
> Coefficient
I'm too old, I didn't find it in the 5th edition, but if Stata uses
it, then it must be ok.
Josef
Skipper Seabold
Well, my thinking was that it is _only_ asymptotically normal due to
CLT. So if you only have a small sample then the approximation is just
bad either way - unless you *know* the errors are normal (right...).
You're gonna need a different hammer anyway.
>
> The debate is a bit like what should be the default ddof for numpy.var.
>
> I don't really care much, and I'm fine with using the normal
> distribution as default everywhere, but I like the option to change
> use_t=True
> (which also has the advantage that I could do a quick MonteCarlo with
> a given design matrix and see which is better.)
You could also do
stats.t.sf([res.tvalues for res in results_list], df)
or use a generator.
But ok, good. I'm going to switch this over and push some discrete
changes shortly.
>
> Using the normal by default, but also be a better signal to users that
> we don't have an exact small sample distribution.
>
> I never checked whether it would make much of a difference except for
> very small samples.
>
> (bootstrap or other approximation might be more important in the small
> sample case anyway.)
>
>
>
> >
> > This leads me to think that pvalues and conf_int in
> > LikelihoodModelResults should all use the normal distribution, and we
> > should overwrite these methods in RegressionResults. Can we agree on
> > this? I need to go back and replace the use of t in some of the other
> > models (it's still incorrectly t in ARMA in addition to where you've
> > pointed out)
>
> I still think the LikelihodModelResults conf_int should base the
> results on the use_t or use_norm if we want to change the name, we can
> have a generic top-level implementation, and it can be changed case by
> case.
So right now this would be use_t = False at the top-level and use_t =
True in linear_model. I see where this *might* be desirable down the
road to avoid any code duplication between conf_int when it takes
method, but I leave this change for when that gets fleshed out. Right
now they're separated to have different notes sections in the docs. It
might be the case that conf_int should be a general function with a
method calling it and giving a distribution under the hood even.
josef...@gmail.com
One thought: for consistency, we should introduce a z_test in analogy
to the t_test method. Currently (or when it is changed to normal), we
get different distribution reported in the pvalues, then in
t_test(np.eye(len(res.params)))
Using the normal distribution directly in t_test will be too much of a
surprise given the name.
What is the normal distribution analog to the f_test for joint linear
restrictions? Wald ?
We need to write the docs if we get the full menagerie of tests.
Josef
Skipper Seabold
<snip>
> One thought: for consistency, we should introduce a z_test in analogy
> to the t_test method. Currently (or when it is changed to normal), we
> get different distribution reported in the pvalues, then in
> t_test(np.eye(len(res.params)))
>
> Using the normal distribution directly in t_test will be too much of a
> surprise given the name.
>
> What is the normal distribution analog to the f_test for joint linear
> restrictions? Wald ?
My impression is that t_test is a Wald test. Stata uses the generic
"test" statement, which uses a chi-2 distribution if the estimator
uses the standard normal and F if t-distribution for joint hypotheses.
Greene 5.6 7th edition, "Nonnormal Disturbances and Large-Sample Tests
[for the linear regression model]" confuses the issue even further.
This is for Wald tests of restrictions on a single coefficient. "The
appropriate critical values only *converge* to those from the standard
normal, and generally *from above*, although we cannot be sure of
this. In the interest of conservatism -- that is, in controlling the
probability of Type I error - one should generally use the critical
values from the t distribution even in the absence of normality."
There is a similar conclusion for the F test. In practice, it should
rarely matter and if you're changing conclusions with a sample of 35
based on a p-value of .058 vs .05, you have bigger problems to worry
about! I think providing t_test, f_test, z_test, and chi2_test should
be fine for completeness, clarity, and transparency. With a formula
framework, we can hide some of this from the user.
results.test("employment = 0 & education = 0") - does a chi2-test
under the hood and reports it as such, or whatever is appropriate (and
hopefully consistent! no pun intended).
Sheesh.
Skipper Seabold
For posterity's sake, there are more details of asymptotic normality
of robust estimators in van der Vaart's Asymptotic Statistics -
Chapter 5, in particular.
Skipper
(from down the rabbit hole)
josef...@gmail.com
nice chapter, thanks for the reference
it has a good list of examples for general M-estimators, not just
robust, and the first time I see the asymptotic normality for the
general case.
Josef