Robust SE & Heteroskedasticity-consistent estimation

An embedded and charset-unspecified text was scrubbed...
Name: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20100510/fce71906/attachment.pl>

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 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

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?
Ah, I didn't realize that the request could have been about the "maxLik" 
package (as opposed to generic maximum likelihood estimation).

Yes, bread() should be easy and estfun() is also easy provided you have 
analytical gradients.

The bread() method can simply return
   vcov(x) * number_of_observations

The estfun() method should return an n x k matrix (n = #obs, k = #par) 
where each row contains the gradient contribution of the corresponding 
observation. Thus
   colSums(estfun(x))
would correspond to
   x$gradient
You typically only get the latter from numerical gradients, but with 
analytical gradients you can easily evaluate it observation-wise.

With these two methods available, you can then call
   sandwich(x)
   sandwich(x, adjust = TRUE)
which computes the basic robust standard errors. In HC terminology, these 
are the HC0 and HC1 estimators.

Calling vcovHC() is more difficult because you have to assure that the 
model is based on a single linear predictor. Also I'm not sure about the 
applicability of HC2-HC4 outside the linear regression model.

For more details, see the "sandwich" vignettes.

Best,
Z
/Arne

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

Thank you all of you,;
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);

n<-NROW(v)
k<-NROW(resMaxlik$estimate)

meat <-function(obj,nrow,npar,adjust=FALSE,...) 
{ psi<-obj$gradient
  k<-npar
  n<-nrow
  rval<-crossprod(as.matrix(psi))/n
  if (adjust) rval <- n/(n-k)*rval
  rval
}

M<-as.vector(meat(resMaxlik,n,k))
B <-as.matrix(ginv(resMaxlik$hessian))

vcovhc=B*M*B

#new standard errors
coeftest(resMaxlik,vcov.=vcovhc)

do you think it's sound good ?

-----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 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
Hi,
I think my "meat"function is wrong.
But the "estfun(resMaxlik)" doesn't work;
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?: RATIARISON Eric 
Envoy??: mardi 11 mai 2010 11:47
??: 'Arne Henningsen'; Achim Zeileis; r-help at r-project.org; Ott-Siim Toomet
Objet?: RE: [R] Robust SE & Heteroskedasticity-consistent estimation

Thank you all of you,;
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);

n<-NROW(v)
k<-NROW(resMaxlik$estimate)

meat <-function(obj,nrow,npar,adjust=FALSE,...) 
{ psi<-obj$gradient
  k<-npar
  n<-nrow
  rval<-crossprod(as.matrix(psi))/n
  if (adjust) rval <- n/(n-k)*rval
  rval
}

M<-as.vector(meat(resMaxlik,n,k))
B <-as.matrix(ginv(resMaxlik$hessian))

vcovhc=B*M*B

#new standard errors
coeftest(resMaxlik,vcov.=vcovhc)

do you think it's sound good ?

-----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 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

Hi,
I think my "meat"function is wrong.
Yes.
But the "estfun(resMaxlik)" doesn't work;
Because you haven't supplied one.

It's as I wrote in my previous mails: You need to write a bread() and an 
estfun() method as specified in my previous mail. Details are in the 
sandwich vignette I already pointed you to. The meat() function does not 
have to be modified.
Z
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?: RATIARISON Eric
Envoy??: mardi 11 mai 2010 11:47
??: 'Arne Henningsen'; Achim Zeileis; r-help at r-project.org; Ott-Siim Toomet
Objet?: RE: [R] Robust SE & Heteroskedasticity-consistent estimation

Thank you all of you,;
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);

n<-NROW(v)
k<-NROW(resMaxlik$estimate)

meat <-function(obj,nrow,npar,adjust=FALSE,...)
{ psi<-obj$gradient
 k<-npar
 n<-nrow
 rval<-crossprod(as.matrix(psi))/n
 if (adjust) rval <- n/(n-k)*rval
 rval
}

M<-as.vector(meat(resMaxlik,n,k))
B <-as.matrix(ginv(resMaxlik$hessian))

vcovhc=B*M*B

#new standard errors
coeftest(resMaxlik,vcov.=vcovhc)

do you think it's sound good ?

-----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

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 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
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

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