[R] Robust SE & Heteroskedasticity-consistent estimation

[Ott-Siim Toomet]

Hi,

the error message tells you that there is no method "estfun" for maxLik
objects.  Which is (unfortunately) true.  I try to look at this stuff
during the weekend (given my son sleeps well ;-).

Actually, maxLik can only implement estfun & friends given the user
supplies them.  AFAIK, you need to supply your likelihood function in
"BHHH form", i.e. without 'sum' in the last row in your example below, and
the gradient in the same way.  If your function only returns a single
likelihood value, there is now way maxLik can "distribute" it to
individual observations.

Best,
Ott



> Hi, here my new version:
> I submit you my case:
> ( Pseudo likehood for exponential family with offset )
>
> loglik <- function(param)  {
>    b<-param
>    m=as.vector(of+z%*%b)
>    ll <- sum(v*m-exp(m))   }
>
> gradlik <- function(param) {
>      b<-param
>      m=as.vector(of+z%*%b)
>      gg<-(v-exp(m))*z      }
>
> hesslik <- function(param) {
>      b<-param
>      m=as.vector(of+z%*%b)
>      hh <- -t(exp(m)*z)%*%z }
>
> resMaxlik <-
> maxLik(loglik,grad=gradlik,hess=hesslik,start=sc,method="nr",tol=1e-4);
>
> resMaxlik$offset<-of
> resMaxlik$x<-z
> resMaxlik$y<-v
>
>
>
> estfun <- function(obj,...)
> {
> m=as.vector(obj$offset+obj$x%*%obj$estimate)
> (obj$y-exp(m))*obj$x
> }
>
> M <- function(obj, adjust = FALSE, ...)
> {
> psi <- estfun(obj)
> k <- NCOL(psi)
> n <- NROW(psi)
> rval <- crossprod(as.matrix(psi))/n
> if(adjust) rval <- n/(n - k) * rval
> rval
> }
>
>
> B <- function(obj, ...)
> { as.matrix(ginv(obj$hessian))
> }
>
>
> So S <-B(resMaxlik)*M(resMaxlik)*B(resMaxlik) directly. It's ok.
>
> But I call sandwich function like this :
> sandwich(resMaxlik,meat.=M,bread.=B)
> It returns a error message : "Erreur dans UseMethod("estfun") :
>   pas de méthode pour 'estfun' applicable pour un objet de classe
> "c('maxLik', 'maxim', 'list')"
>
>
>
>
>
> -----Message d'origine-----
> De : Arne Henningsen [mailto:arne.henningsen at googlemail.com]
> Envoyé : mardi 11 mai 2010 08:25
> À : Achim Zeileis; RATIARISON Eric; r-help at r-project.org; Ott-Siim Toomet
> Objet : Re: [R] Robust SE & Heteroskedasticity-consistent estimation
>
> On 11 May 2010 00:52, Achim Zeileis <Achim.Zeileis at uibk.ac.at> wrote:
>> On Mon, 10 May 2010, RATIARISON Eric wrote:
>>> I'm using maxlik with functions specified (L, his gradient & hessian).
>>>
>>> Now I would like determine some robust standard errors of my
>>> estimators.
>>>
>>> So I 'm try to use vcovHC, or hccm or robcov for example
>>>
>>> but in use one of them with my result of maxlik, I've a the following
>>> error message :
>>>
>>> Erreur dans terms.default(object) : no terms component
>>>
>>> Is there some attributes to give to maxlik objet for "fitting" the call
>>> of vcovHC?
>>
>> This is discussed in
>>  vignette("sandwich-OOP", package = "sandwich")
>> one of the vignettes accompanying the "sandwich" package that provides
>> the
>> vcovHC() function. At the very least, you need an estfun() method which
>> extracts the gradient contributions per observation. Then you need a
>> bread()
>> function, typically based on the observed Hessian. Then you can compute
>> the
>> basic sandwich() estimators.
>
> Is it possible to implement the estfun() method and the bread()
> function in the maxLik package so that vcovHC() can be easily used by
> all users of the maxLik package? If yes: should we (Eric, Achim, Ott,
> Arne) implement this feature together?
>
> /Arne
>
> --
> Arne Henningsen
> http://www.arne-henningsen.name
>