[R] poisson regression with robust error variance ('eyestudy')

[Frank E Harrell Jr]

Lutz Ph. Breitling wrote:
> Thank you very much for your comments!
> 
> however, i still do not get it right.
> robcov() accepts fit objects like lrm or ols objects as arguments,
> but obviously not the glmD objects (or at least not as simple as that).
> below some code to demonstrate.
> what am i still doing wrong?
> 
> thx for your efforts-
> lutz
> 
> 
> id<-1:500
> outcome<-sample(c(0,1), 500, replace=T, prob=c(.6, .4))
> exposed<-sample(c(0,1), 500, replace=T, prob=c(.5, .5))
> my.data<-data.frame(id=id, ou=outcome, ex=exposed)
> 
> model1<-glmD(ou~ex, my.data, family=poisson(link="log"))
> robcov(model1)
> 
> Error in match.arg(type) : ARG should be one of deviance, pearson, 
> working, response, partial

Sorry I didn't think of that sooner.  robcov needs the residuals method 
for the fitter to allow a type="score" or type="hscore" (for Efron's 
method) argument.  Until someone adds score residuals to residuals.glm 
robcov will not work for you.  residuals.lrm and residuals.coxph are 
examples where score residuals are computed.  You can get robust 
variance-covariance estimates with the bootstrap using bootcov for glmD 
fits.  Oddly in your example I am finding that the bootstrap variances 
are lower than the information-matrix-based ones.

Frank Harrell

> 
> At 17:25 02.06.2004, Frank E Harrell Jr wrote:
> 
>> Thomas Lumley wrote:
>>
>>> On Wed, 2 Jun 2004, Lutz Ph. Breitling wrote:
>>>
>>>> Dear all,
>>>>
>>>> i am trying to redo the 'eyestudy' analysis presented on the site
>>>> http://www.ats.ucla.edu/stat/stata/faq/relative_risk.htm
>>>> with R (1.9.0), with special interest in the section on "relative risk
>>>> estimation by poisson regression with robust error variance".
>>>>
>>>> so i guess rlm is the function to use. but what is its equivalent to 
>>>> the
>>>> glm's argument "family" to indicate 'poisson'? or am i somehow totally
>>>> wrong and this is not applicable here?
>>>
>>>
>>> No, no.  You want glm() and then a function to compute the robust
>>> covariance matrix (there's robcov() in the Hmisc package), or use gee()
>>> from the "gee" package or geese() from "geepack" with independence 
>>> working
>>> correlation.
>>
>>
>> Slight correction: robcov in the Design package, can easily be used 
>> with Design's glmD function.  -Frank
>>
>>> These are not outlier-resistant estimates of the regression 
>>> coefficients,
>>> they are model-agnostic estimates of the standard errors.
>>> Stata is unusual in providing these covariance matrix estimates for just
>>> about every regression estimator.  I think R should consider doing
>>> something similar, but haven't got around to it.
>>>         -thomas
>>>
>>>
>>>> thx a lot-
>>>> lutz
>>>>
>>>>
>>>> =============================
>>>> Lutz Ph. Breitling, CMd
>>>
>>>
>>> Thomas Lumley                   Assoc. Professor, Biostatistics
>>> tlumley at u.washington.edu        University of Washington, Seattle
> 
> =============================
> Lutz Ph. Breitling
> Unité des Recherches Médicale
> Hôpital Albert Schweitzer
> B.P. 118 Lambaréné (GABON)
> 
> 

-- 
Frank E Harrell Jr   Professor and Chair           School of Medicine
                      Department of Biostatistics   Vanderbilt University