R help - May 2010 - Robust SE & Heteroskedasticity-consistent estimation

Hi,

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?

 

Thank you

 


	[[alternative HTML version deleted]]

On Mon, 10 May 2010, RATIARISON Eric wrote:
> Hi,
>
> 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.

Best,
Z
>
>
> Thank you
>
>
>
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

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