Error in logistic regression with "cloglog" link

[Thomas Haslwanter]

In trying to reproduce the code results from Dobson & Barnett "An
Introduction to GLMs" in Python, I have run into an inconsistency between
the results from R and from Python. While the results for logistic
regression with statsmodels match the R-results for the logit and probit
link functions, the results for the cloglog link are inconsistent.

I have looked at the Python code in statsmodels, and it seems correct to me,
so I am a bit dumbfounded.
Also, the statsmodels link only works for "cloglog", but crashes for
"CLogLog".


Python Code:
------------
'''Solution to Exercise 7.3.1 (Beetle Mortality) in Dobson & Bartlett:
"Introduction to GLMs"

author: thomas haslwanter
date:   may 2013
ver:    1.0

'''

# Standard libraries
import numpy as np
import pandas as pd
import statsmodels.api as sm
   
if __name__ == '__main__':
    data = np.array([[1.6907, 59, 6],
        [1.7242, 60, 13],
        [1.7552, 62, 18],
        [1.7842, 56, 28],
        [1.8113, 63, 52],
        [1.8369, 59, 53],
        [1.861, 62, 61 ],
        [1.8839, 60, 60]])
   
    df = pd.DataFrame(data)
    df.columns = ['x', 'n', 'y']
   
    # fit the model
    df['tested'] = df['n']
    df['killed'] = df['y']
    df['survived'] = df['tested'] - df['killed']
   
    model_logit = sm.GLM.from_formula('survived + killed ~ x', data=df, family=sm.families.Binomial(sm.families.links.logit)).fit()
    print 'Logit:'
    print df['n']*(1-model_logit.fittedvalues)
   
    model_cll = sm.GLM.from_formula('survived + killed ~ x', data=df, family=sm.families.Binomial(sm.families.links.cloglog)).fit()
    print 'CLogLog:'
    print df['n']*(1-model_cll.fittedvalues)

Python Results:
---------------
Logit:
0     3.457461
1     9.841672
2    22.451378
3    33.897635
4    50.095822
5    53.290913
6    59.222159
7    58.742961
CLogLog:
0     2.365465
1    12.557569
2    27.787944
3    36.431779
4    49.561865
5    51.380562
6    57.100959
7    57.057706

R code:
-------
y=c(6,13,18,28,52,53,61,60)
n=c(59,60,62,56,63,59,62,60)
x=c(1.6907,1.7242,1.7552,1.7842,1.8113,1.8369,1.8610,1.8839)
n_y=n-y
beetle.mat=cbind(y,n_y)
res.glm1=glm(beetle.mat~x, family=binomial(link="logit"))
res.glm3=glm(beetle.mat~x, family=binomial(link="cloglog"))
fitted.values(res.glm1)*n
fitted.values(res.glm3)*n

R results:
----------
> fitted.values(res.glm1)*n
        1         2         3         4         5         6         7         8
 3.457461  9.841672 22.451378 33.897635 50.095822 53.290913 59.222159 58.742961
> fitted.values(res.glm3)*n
       1        2        3        4        5        6        7        8
 5.58945 11.28068 20.95422 30.36944 47.77642 54.14273 61.11331 59.94723

[josef...@gmail.com]

Are params for cloglog also different from R?

Josef

[josef...@gmail.com]

It looks like cloglog link has p and 1-p reversed

model_cll2 = sm.GLM.from_formula('killed + survived ~ x', data=df,
family=sm.families.Binomial(sm.families.links.cloglog)).fit()
print 'CLogLog:'
print df['n']*(model_cll2.fittedvalues)

code comment says

#TODO: CLogLog is untested
class CLogLog(Logit):

and CLogLog is not registered as valid link, it seems.

Please open an issue?


Josef

[josef...@gmail.com]

the problem might be that we need a loglog link
cloglog looks correct for complementary log log

http://www.econ.uiuc.edu/~roger/research/links/links.pdf

VGAM has loglog and cloglog
http://rss.acs.unt.edu/Rdoc/library/VGAM/html/cloglog.html

My guess from a quick browsing is that it depends on whether p or 1-p
is very small

Josef


>
>
> Josef

[josef...@gmail.com]

http://r.789695.n4.nabble.com/glm-binomial-loglog-NOT-cloglog-link-td849370.html

>
> Josef
>
>
>>
>>
>> Josef

[josef...@gmail.com]

And here is something why R's cloglog is actually loglog

http://tolstoy.newcastle.edu.au/R/help/99a/0237.htm
http://tolstoy.newcastle.edu.au/R/help/99a/0238.htm

Josef

>
>>
>> Josef
>>
>>
>>>
>>>
>>> Josef

[josef...@gmail.com]

another fun one Stata versus SPSS

http://statalist.blogspot.ca/2006/01/st-re-cloglog-loglog-etc.html

Josef
<back to other things>

>
> Josef
>
>>
>>>
>>> Josef
>>>
>>>
>>>>
>>>>
>>>> Josef

[Skipper Seabold]

On Sat, May 11, 2013 at 10:45 AM, <josef...@gmail.com> wrote:
> another fun one Stata versus SPSS
>
> http://statalist.blogspot.ca/2006/01/st-re-cloglog-loglog-etc.html
>

Hardin and Hilbe are a big reason I chose to test against Stata vs. R
for GLMs (except cloglog...). They are tough to beat.

Skipper