Disagreement of confidence intervals between OLS and GLM

[Thomas Haslwanter]

Yet another comment - I hope I am not becoming a pain in the neck ;)

When I calculate model parameters with OLS and with GLM, I get the same parameters - but I obtain different values for the confidence intervals (CIs)!
When I do the CIs by hand, the values I get are equivalent to those reported by OLS.
My question: I think that the CIs should be the same - should they not?
Also: after fitting an OLS model, how can I obtain the standard errors displayed by "model.summary()"?

I have the code I have written so far (Python/statsmodels solutions to the Dobson book on GLMs) on github
   g...@github.com:thomas-haslwanter/dobson.git
(or the link https://github.com/thomas-haslwanter/dobson)

and you can see the differences simply by running "dobson.py". The lines producing the same parameters, but different CIs are

model = sm.OLS.from_formula('carbohydrate ~ age + weight + protein', data=df).fit()
print model.summary()
glm = sm.GLM.from_formula('carbohydrate ~ age + weight + protein', family=sm.families.Gaussian(), data=df).fit()
print glm.summary()

and the differing outputs are

                            OLS Regression Results                           
==============================================================================
Dep. Variable:           carbohydrate   R-squared:                       0.481
Model:                            OLS   Adj. R-squared:                  0.383
Method:                 Least Squares   F-statistic:                     4.934
Date:                Sun, 12 May 2013   Prob (F-statistic):             0.0130
Time:                        18:23:09   Log-Likelihood:                -61.837
No. Observations:                  20   AIC:                             131.7
Df Residuals:                      16   BIC:                             135.7
Df Model:                           3                                        
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept     36.9601     13.071      2.828      0.012         9.250    64.670
age           -0.1137      0.109     -1.040      0.314        -0.345     0.118
weight        -0.2280      0.083     -2.738      0.015        -0.405    -0.051
protein        1.9577      0.635      3.084      0.007         0.612     3.304
==============================================================================
Omnibus:                        2.146   Durbin-Watson:                   1.875
Prob(Omnibus):                  0.342   Jarque-Bera (JB):                1.425
Skew:                          -0.417   Prob(JB):                        0.490
Kurtosis:                       1.993   Cond. No.                     1.20e+03
==============================================================================

Same model, calculated with GLM
                 Generalized Linear Model Regression Results                 
==============================================================================
Dep. Variable:           carbohydrate   No. Observations:                   20
Model:                            GLM   Df Residuals:                       16
Model Family:                Gaussian   Df Model:                            3
Link Function:               identity   Scale:                   35.4789285829
Method:                          IRLS   Log-Likelihood:                -61.837
Date:                Sun, 12 May 2013   Deviance:                       567.66
Time:                        18:23:09   Pearson chi2:                     568.
No. Iterations:                     3                                        
==============================================================================
                 coef    std err          t      P>|t|      [95.0% Conf. Int.]
------------------------------------------------------------------------------
Intercept     36.9601     13.071      2.828      0.012        11.341    62.579
age           -0.1137      0.109     -1.040      0.314        -0.328     0.101
weight        -0.2280      0.083     -2.738      0.015        -0.391    -0.065
protein        1.9577      0.635      3.084      0.007         0.713     3.202
==============================================================================

[josef...@gmail.com]

On Sun, May 12, 2013 at 12:25 PM, Thomas Haslwanter
<thomas.h...@gmail.com> wrote:
> Yet another comment - I hope I am not becoming a pain in the neck ;)
>
> When I calculate model parameters with OLS and with GLM, I get the same
> parameters - but I obtain different values for the confidence intervals
> (CIs)!

quick guess: one uses the t-distribution (OLS) the other uses the
normal distribution (GLM).
In small samples the difference can be large.

(I can look at the other questions in a bit)

Josef

[Skipper Seabold]

On Sun, May 12, 2013 at 12:36 PM, <josef...@gmail.com> wrote:
> On Sun, May 12, 2013 at 12:25 PM, Thomas Haslwanter
> <thomas.h...@gmail.com> wrote:
>> Yet another comment - I hope I am not becoming a pain in the neck ;)
>>
>> When I calculate model parameters with OLS and with GLM, I get the same
>> parameters - but I obtain different values for the confidence intervals
>> (CIs)!
>
> quick guess: one uses the t-distribution (OLS) the other uses the
> normal distribution (GLM).
> In small samples the difference can be large.
>

Yes, OLS uses the t and GLM uses the normal IIRC.

Skipper

[Nathaniel Smith]

Which means: OLS uses a method that gives exact results, but only works in the special case where all the usual OLS criteria apply - iid Gaussian noise etc. GLM instead uses an approximate method which is correct asymptotically but may be off for small samples; the tradeoff you get in return is that this method works the same way for all GLM models, including those with non-Gaussian error terms and non-trivial link functions. So that's why they're different.

-n

[Thomas Haslwanter]

Thanks, Skipper - I am impressed how quickly you figure such things out!

[josef...@gmail.com]

On Sun, May 12, 2013 at 12:25 PM, Thomas Haslwanter
<thomas.h...@gmail.com> wrote:

> Yet another comment - I hope I am not becoming a pain in the neck ;)

I think it's very useful what you are doing.
And, if you find pieces that are missing, then you could open issues
for them, so we know which gaps need filling.

>
> When I calculate model parameters with OLS and with GLM, I get the same
> parameters - but I obtain different values for the confidence intervals
> (CIs)!
> When I do the CIs by hand, the values I get are equivalent to those reported
> by OLS.
> My question: I think that the CIs should be the same - should they not?
> Also: after fitting an OLS model, how can I obtain the standard errors
> displayed by "model.summary()"?

bse, tvalues and conf_int() should give you the same numbers as the summary
they should be basic attributes across models.
tvalues could be zvalues if we use normal distribution, but we decided
to use the same name for both for consistency.

for OLS (and example I had open)
>>> result.bse
array([ 0.14488718, 0.16562277, 0.31758337])
>>> result.tvalues
array([ 5.74811132, 12.32127366, 5.88324762])
>>> result.conf_int()

I think it's a mistake that it shows t and P>|t| instead of z and P>|z|

We changed whether we use t or normal distribution in the different
models a while ago after some discussion.
I don't know if there is a way for the user to switch.
I need to check this.

Josef

[josef...@gmail.com]

I think you found a bug in GLMResults.

pvalues still use the t-distribution, conf_int use the normal distribution
Looks like some incomplete refactoring that escaped us.

There is currently no way to overwrite which distribution is used.

Josef


>
> Josef

[josef...@gmail.com]

It's not always very clear whether t or normal is better.
With large numbers they are the same.
In small samples some textbooks or articles report that t and F
distributions work better than normal and chisquare distributions even
if we get the (asymptotic) theory only for the latter.

a lengthy discussion
https://groups.google.com/d/msg/pystatsmodels/5Y7pUYR8Ib4/ObrbbWUOw0YJ
maybe not the only one we had.
But I don't find a paper trail on github. Maybe it was changed before
we started to use pull request for every change.

Josef

>
> -n