[R] Understanding R code

[David Winsemius]

On Aug 20, 2009, at 4:31 PM, kfcnhl wrote:

>
> What is
> 1. par.ests <- optimfit$par

Just a guess... the parameter estimates from an optimization procedure?

> 2. fisher <- hessb(negloglik, par.ests, maxvalue=maxima);

Probably the Fisher Information matrix.

> 3. varcov <- solve(fisher);

I seem to remember being told that the inverse of the Fisher  
Information matrix gives you the variance covariance matrix.

> 4. par.ses <- sqrt(diag(varcov));

And the square roots of the diagonal elements of the variance- 
covariance matrix would be the standard errors of the parameter  
estimates.

Is the pop quiz over?


>
> Thanks a lot,
>
> fit.GEV <- function(maxima)
> {
> sigma0 <- sqrt((6. * var(maxima))/pi)
> mu0 <- mean(maxima) - 0.57722 * sigma0
> xi0 <- 0.1
> theta <- c(xi0, mu0, sigma0)
>
> #10/5/2007: removed assign() for maxima.nl
> #10/5/2007: passed additional parameter to avoid using assign()
> negloglik <- function(theta, maxvalue)
> {
> -sum(dGEV(maxvalue,theta[1],theta[2],abs(theta[3]),logvalue=TRUE));
> }
>
> #10/5/2007: passed additional 4th parameter which will be passed to
> negloglik
> optimfit <- optim(theta, fn=negloglik, gr=NULL, maxvalue=maxima,
> method="BFGS");
> par.ests <- optimfit$par;
> if(optimfit$convergence == 0) #if code is 0 we are OK;
> converged <- TRUE
> else
> converged <- FALSE;
> #replace 3rd parameter estimate with its absolute value:
> par.ests[3] <- abs(par.ests[3]);
> #10/5/2007: added final parameter which must be passed to negloglik
> fisher <- hessb(negloglik, par.ests, maxvalue=maxima);
> varcov <- solve(fisher);
>
> par.ses <- sqrt(diag(varcov));
> out <- list(par.ests = par.ests, par.ses = par.ses, varcov = varcov,
> converged =
> #10/5/2007: passed additional 2nd parameter to -negloglik()
> converged, llmax = -negloglik(par.ests,maxima));
> #06/30/2008: the order of the names is incorrect; note that the  
> parameter
> vector theta
> #(passed to the optimization function) has the order c(xi0, mu0,  
> sigma0);
> the order of
> #the return parameters should be identical;
> #hence we need to make the following correction to maintain the same  
> order
> #names(out$par.ests) <- c("xi", "sigma", "mu");
> #names(out$par.ses) <- c("xi", "sigma", "mu");
> names(out$par.ests) <- c("xi", "mu", "sigma");
> names(out$par.ses) <- c("xi", "mu", "sigma");
> out;
> }
> -- 
> View this message in context: http://www.nabble.com/Understanding-R-code-tp25069337p25069337.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> 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.

David Winsemius, MD
Heritage Laboratories
West Hartford, CT