Contour2

Here is an example of using the function “contour” from the base graphics package in R to draw a contour plot. This improves and extends a previous version.

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)
Nobs <- 1000
Data <- rgamma(Nobs, shape = 6, rate = 0.5)
library(MASS)
truehist(Data, xlim = c(0, 35))
xpts <- seq(from = 0, to = 35, length = 1000)
lines(xpts, dgamma(xpts, shape = 6, rate = 0.5))

fit <- fitdistr(Data, "gamma")

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

fit

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

NsPts <- 100
NrPts <- 100

sGridPts <- seq(from = 5, to = 7, length = NsPts)
rGridPts <- seq(from = 0.4, to = 0.6, length = NrPts)

triples <- expand.grid(list(sGridPts = sGridPts, rGridPts = rGridPts, Data = Data))
       # All combinations of values of shape and rate parameters from 
       # the two grids, and an observation

head(triples)

##   sGridPts rGridPts  Data
## 1    5.000      0.4 6.366
## 2    5.020      0.4 6.366
## 3    5.040      0.4 6.366
## 4    5.061      0.4 6.366
## 5    5.081      0.4 6.366
## 6    5.101      0.4 6.366

str(triples)

## 'data.frame':    10000000 obs. of  3 variables:
##  $ sGridPts: num  5 5.02 5.04 5.06 5.08 ...
##  $ rGridPts: num  0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 ...
##  $ Data    : num  6.37 6.37 6.37 6.37 6.37 ...
##  - attr(*, "out.attrs")=List of 2
##   ..$ dim     : Named int  100 100 1000
##   .. ..- attr(*, "names")= chr  "sGridPts" "rGridPts" "Data"
##   ..$ dimnames:List of 3
##   .. ..$ sGridPts: chr  "sGridPts=5.000" "sGridPts=5.020" "sGridPts=5.040" "sGridPts=5.061" ...
##   .. ..$ rGridPts: chr  "rGridPts=0.4000" "rGridPts=0.4020" "rGridPts=0.4040" "rGridPts=0.4061" ...
##   .. ..$ Data    : chr  "Data= 6.366" "Data=10.732" "Data= 4.837" "Data=15.767" ...

LogLikOnGrid <- apply(
    X = array(dgamma(triples$Data, 
            shape = triples$sGridPts, 
            rate = triples$rGridPts, 
            log = TRUE), 
                dim = c(NsPts, NrPts, Nobs), 
                dimnames = attr(triples, "out.attrs")$dimnames),
    MARGIN = c("sGridPts", "rGridPts"),
    FUN = sum)



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


MaxLogLik <- max(LogLikOnGrid)
LogLikOnGrid <- matrix(LogLikOnGrid, NsPts, NrPts) - max(LogLikOnGrid)




contour(x = sGridPts, y = rGridPts, z = LogLikOnGrid)

contour(x = sGridPts, y = rGridPts, z = LogLikOnGrid, 
        levels = c( -(1:8)*5), col = "blue", add = TRUE)

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

## Just to show the wire-frame plot

persp(x = sGridPts, y = rGridPts, z = LogLikOnGrid)

## How to do this using the "lattice" package


pairs <- expand.grid(list(sGridPts = sGridPts, rGridPts = rGridPts))

newTriples <- cbind(pairs, LogLikOnGrid = as.vector(LogLikOnGrid))

head(newTriples)

##   sGridPts rGridPts LogLikOnGrid
## 1    5.000      0.4       -14.49
## 2    5.020      0.4       -15.11
## 3    5.040      0.4       -15.83
## 4    5.061      0.4       -16.63
## 5    5.081      0.4       -17.53
## 6    5.101      0.4       -18.51

library(lattice)

library(latticeExtra)

## Loading required package: RColorBrewer

Part1 <- contourplot(LogLikOnGrid ~ sGridPts * rGridPts, data = newTriples, 
    cuts = 10)

Part2 <- contourplot(LogLikOnGrid ~ sGridPts * rGridPts, data = newTriples, 
        at = c( -(1:8)*5), col = "blue")

Part3 <- contourplot(LogLikOnGrid ~ sGridPts * rGridPts, data = newTriples, 
        at = - qchisq(0.95, 2)/2, col = "red")
        
print(Part1 + Part2 + Part3)