[ADMB Users] extract Hessian from within R

[Ben Bolker]

  Since this was cc'd to me:  I would recommend using read_pars() from
the R2admb package, which has all of this output-reading stuff (and
more) built in.  If you don't want to do that, you can go to

https://github.com/bbolker/R2admb/blob/master/R2admb/R/read-funs.r

and find the code for reading the Hessian file:

 if(all(is.na(vcov)) & file.exists("admodel.hes"))
    {
        filen <- file("admodel.hes", "rb")
        nopar <- readBin(filen, what = "integer", n = 1)
        hes <- readBin(filen, what = "double", n = nopar * nopar)
        hes <- matrix(hes, byrow = TRUE, ncol = nopar)
        colnames(hes) <- rownames(hes) <- sdparnames[seq(ncol(vcov))]
        close(filen)
    }

  cheers
    Ben Bolker


On 13-11-04 04:40 AM, Allan Clark wrote:
> Hello all
> 
> I've coded up a simple model where I estimate various probabilities.
> These are restricted to (0,1) and decided to rather transform these
> variables so that they are not bounded anymore. How does one extract the
> hessian of the parameters from within R?
> 
> i've found functions on the following cite useful:
> http://qfc.fw.msu.edu/userfun.asp#r3 . But are there other wasys of
> extracting the Hessian as well?
> 
> (Once one has the hessian in terms of the transformed parameters, it is
> eay to calculate the hessian in terms of the initial variables.)
> 
> regards
> 
> 
> 
> 
> 
> read.admbFit()
> 
> |
> ## Function to read AD Model Builder fit from its par and cor files.
> ## Use by:
> ## simple.fit <- read.admbFit('c:/admb/examples/simple')
> ## Then the object 'simple.fit' is a list containing sub-objects
> # 'names', 'est', 'std', 'cor', and 'cov' for all model
> # parameters and sdreport quantities.
> #
> read.admbFit<-function(file){
>   ret<-list()
>   # read par file
>   parfile<-as.numeric(scan(paste(file,'.par', sep=''),
>     what='', n=16, quiet=TRUE)[c(6,11,16)])
>   ret$nopar<-as.integer(parfile[1])
>   ret$nloglike<-parfile[2] #objective function value
>   ret$maxgrad<-parfile[3]
>   
>   # read cor file
>   file<-paste(file,'.cor', sep='')
>   lin<-readLines(file)
>   # total parameter including sdreport variables
>   ret$totPar<-length(lin)-2
>   #log of the determinant of the hessian
>   ret$logDetHess<-as.numeric(strsplit(lin[1], '=')[[1]][2])
>   sublin<-lapply(strsplit(lin[1:ret$totPar+2], ' '),function(x)x[x!=''])
>   ret$names<-unlist(lapply(sublin,function(x)x[2]))
>   ret$est<-as.numeric(unlist(lapply(sublin,function(x)x[3])))
>   ret$std<-as.numeric(unlist(lapply(sublin,function(x)x[4])))
>   ret$cor<-matrix(NA, ret$totPar, ret$totPar)
>   corvec<-unlist(sapply(1:length(sublin), function(i)sublin[[i]][5:(4+i)]))
>   ret$cor[upper.tri(ret$cor, diag=TRUE)]<-as.numeric(corvec)
>   ret$cor[lower.tri(ret$cor)] <- t(ret$cor)[lower.tri(ret$cor)]
>   # covariance matrix
>   ret$cov<-ret$cor*(ret$std %o% ret$std)
>   return(ret)
> }
> |
> 
>