Contour

Here is an example of using the function “contour” from the base graphics package in R to draw a contour plot.

I will generate a random sample of size 1000 from a gamma distribution with certain shape and rate parameters. I’ll estimate the parameters using maximum likelihood, there is a very convenient function for this from the package “MASS” called “fitdistr”.

I’ll also draw the contours of the log-likelihood surface and add the boundary of an approximately 95% confidence region for the two parameters determined using the Wilks “minus twice log likelihood” chi-square approximation.

As you’ll see, the visualisation of the log likelihood function adds hugely to the information which we get from “fitdistr”. And the old-fashioned contour plot seems in this case a much better tool than the more “sexy” wire-frame plot.

One could go on to rotate the surface and even make a movie of it rotating using the package “rgl” but it seems to me that the contour plot is the way to go.

set.seed(110951)
x <- rgamma(1000, 6, 0.5)
library(MASS)
truehist(x, xlim = c(0, 35))
xpts <- seq(from = 0, to = 35, length = 1000)
lines(xpts, dgamma(xpts, 6, 0.5))

fit <- fitdistr(x, "gamma")

## Warning: NaNs produced
## Warning: NaNs produced

fit

##     shape     rate  
##   6.04952   0.50731 
##  (0.26342) (0.02303)

NxPts <- 100
NyPts <- 100

xGridPts <- seq(from = 5, to = 7, length = NxPts)
yGridPts <- seq(from = 0.4, to = 0.6, length = NyPts)

xypairs <- expand.grid(xGridPts, yGridPts) 
       # coordinates of all gridpoints, rows of an N x 2 matrix
       # N = NxPts * NyPts

head(xypairs)

##    Var1 Var2
## 1 5.000  0.4
## 2 5.020  0.4
## 3 5.040  0.4
## 4 5.061  0.4
## 5 5.081  0.4
## 6 5.101  0.4

str(xypairs)

## 'data.frame':    10000 obs. of  2 variables:
##  $ Var1: num  5 5.02 5.04 5.06 5.08 ...
##  $ Var2: num  0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 ...
##  - attr(*, "out.attrs")=List of 2
##   ..$ dim     : int  100 100
##   ..$ dimnames:List of 2
##   .. ..$ Var1: chr  "Var1=5.000" "Var1=5.020" "Var1=5.040" "Var1=5.061" ...
##   .. ..$ Var2: chr  "Var2=0.4000" "Var2=0.4020" "Var2=0.4040" "Var2=0.4061" ...

LogLikFun <- function(shape, rate) sum(dgamma(x, shape = shape, rate = rate, log = TRUE))

LogLikOnGrid <- mapply(LogLikFun, xypairs[ , 1], xypairs[ ,2]) # another long vector

## Can you understand the error message when you do
##    LogLikOnGrid <- outer(xGridPts, yGridPts, LogLikFun)   ?
## Why doesn't "outer" work here??
## Can you come up with something better than my solution?
## You'll also notice that "mapply" is rather slow ...
## There is an explicit "for loop" going on inside of "mapply".
## There ought to be a "vectorised" way to do this.

dim(LogLikOnGrid) <- c(NxPts, NyPts) # convert to matrix


## I'll subtract the maximum from the log likelihood function.
## We are only interested in *differences*, not in its absolute
## value. Now I know the maximum of the surface is at height z = 0.
## Also I reshape as a matrix.

MaxLogLik <- max(LogLikOnGrid)
LogLikOnGrid <- matrix(LogLikOnGrid, NxPts, NyPts) - max(LogLikOnGrid)




contour(x = xGridPts, y = yGridPts, z = LogLikOnGrid)

contour(x = xGridPts, y = yGridPts, z = LogLikOnGrid, 
        levels = c( (0:-8)*10), col = "blue", add = TRUE)

contour(x = xGridPts, y = yGridPts, z = LogLikOnGrid, 
        levels = - qchisq(0.95, 2)/2, col = "red", add = TRUE)

persp(x = xGridPts, y = yGridPts, z = LogLikOnGrid)