[R] Robust SE & Heteroskedasticity-consistent estimation

RATIARISON Eric ERATIARISON at monceauassurances.com
Tue May 18 07:06:46 CEST 2010 

Thanks all of you.


-----Message d'origine-----
De : Arne Henningsen [mailto:arne.henningsen at googlemail.com] 
Envoyé : mardi 18 mai 2010 06:00
À : RATIARISON Eric
Cc : r-help at r-project.org; Ott-Siim Toomet; Achim Zeileis
Objet : Re: [R] Robust SE & Heteroskedasticity-consistent estimation

On 17 May 2010 22:55, RATIARISON Eric <ERATIARISON at monceauassurances.com> wrote:
> It's ok Arne, i've build the MS windows binary.
>
> And the result is ok (sandwich function runs)with this next modifications in my user specified functions
> But with the following functions (individual one):
>
>
> sc= as.vector( c(log(mean(v)),rep(0,dim(z)[2]-1)))
>
> loglik <- function(param)  {
>   b<-param
>   m=as.vector(of+z%*%b)
>   ll <- (v*m-exp(m))   } # no sum
>
> gradlik <- function(param) {
>     b<-param
>     m=as.vector(of+z%*%b)
>     gg <-(v-exp(m))*z      } #here I've replaced %*% by *
>
> 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")
>
> however, I've a question again about bread function. Does it use the  result of hessian defined by user?

Yes.

> When I do all.equal(vcov(resMaxlik),inv(resMaxlik$hessian))
> I've this message : [1] "Mean relative difference: 2" what does it means?

The covariance matrix is the inverse of the *negative* Hessian. So
   all.equal(vcov(resMaxlik),solve(-hessian(resMaxlik)))
should return TRUE

> Lastly, in the last "official" manual "Maxlik.pdf", the sentence "However, no attempt is made to correct the resulting variance-covariance matrix"

Do you mean this paragraph?

     'maxLik' "support" constrained optimization in the sense that
     constraints are passed further to the underlying optimization
     routines, and suitable default method is selected.  However, no
     attempt is made to correct the resulting variance-covariance
     matrix.  Hence, the inference may be wrong.  A corresponding
     warning is issued by the summary method.

This means that the constraints are not taken into account when maxLik
calculates the covariance matrix of the estimates -- the constraints
are simply ignored.

> Should be corrected with possible use of sandwich should'nt it?

I don't think that just using the sandwich method is sufficient to
calculate consistent covariance matrices, when there are (equality or
inequality) constraints on the parameters.

/Arne

> -----Message d'origine-----
> De : Arne Henningsen [mailto:arne.henningsen at googlemail.com]
> Envoyé : lundi 17 mai 2010 17:58
> À : RATIARISON Eric
> Objet : Re: [R] Robust SE & Heteroskedasticity-consistent estimation
>
> On 17 May 2010 17:35, RATIARISON Eric <ERATIARISON at monceauassurances.com> wrote:
>> Thank you Arne,
>> I'm working on PC computer at my work with firewall.
>> How can I reach the new version for testing (in comparing with the gauss procedure named CML) ?
>> I don't understand the link joined.
>
> I have attached the package. You can easily install it if you (also) have a Unix-like OS. If you have MS-Windows, I could build a MS-Windows binary package on http://win-builder.r-project.org/ tomorrow (or you could change the maintainer field and do it yourself).
>
> /Arne
>
>>
>>
>> -----Message d'origine-----
>> De : Arne Henningsen [mailto:arne.henningsen at googlemail.com]
>> Envoyé : lundi 17 mai 2010 17:12
>> À : RATIARISON Eric
>> Cc : Achim Zeileis; r-help at r-project.org; Ott-Siim Toomet Objet : Re:
>> [R] Robust SE & Heteroskedasticity-consistent estimation
>>
>> On 13 May 2010 00:41, RATIARISON Eric <ERATIARISON at monceauassurances.com> wrote:
>>> 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')"
>>
>> I have added an estfun() method and a bread() method for objects of
>> class "maxLik" so that you can use "sandwich" now.
>>
>> Please note:
>>
>> a) The methods can be applied only if maxLik() was called with
>> argument "grad" equal to a gradient function or (if no gradient
>> function is specified) argument "logLik" equal to a log-likelihood
>> function that return the gradients or log-likelihood values,
>> respectively, for *each observation* (i.e. in the BHHH form).
>>
>> b) The package with the new features is not yet available on CRAN but
>> on R-Forge: http://r-forge.r-project.org/scm/?group_id=254
>>
>> Please let me know if you experience any problems.
>>
>> /Arne
>>
>>
>>> -----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
>>>
>>
>>
>>
>> --
>> Arne Henningsen
>> http://www.arne-henningsen.name
>>
>
>
>
> --
> Arne Henningsen
> http://www.arne-henningsen.name
>



-- 
Arne Henningsen
http://www.arne-henningsen.name

More information about the R-help mailing list