%\VignetteIndexEntry{Auckland child mortality 1977-1985 data set}
%\VignetteDepends{}
%\VignetteKeywords{spatial}
%\VignettePackage{spdep}
\documentclass[a4paper,10pt]{article} 
\usepackage{Sweave}
\usepackage{times}
\usepackage{mathptm}
\usepackage{hyperref}
\setkeys{Gin}{width=0.95\textwidth}
\newcommand{\strong}[1]{{\normalfont\fontseries{b}\selectfont #1}}
\let\pkg=\strong
\RequirePackage{alltt}
\newenvironment{example}{\begin{alltt}}{\end{alltt}}
\newenvironment{smallexample}{\begin{alltt}\small}{\end{alltt}}
\newcommand{\code}[1]{\texttt{\small #1}}
\def\RR{\textsf{R}\/}
\def\SP{\texttt{S-PLUS}\/}
\def\SS{\texttt{S}\/}

\title{Auckland child mortality 1977-1985 data set} 
\author{Roger Bivand} 

\begin{document} 

\maketitle 

<<echo=FALSE>>=
owidth <- getOption("width")
options("width"=80)
ow <- getOption("warn")
options("warn"=-1)
.PngNo <- 0
@

<<label=afig,echo=FALSE,eval=FALSE>>= 
.PngNo <- .PngNo + 1; file <- paste("Fig-bitmap-", .PngNo, ".pdf", sep="")
pdf(file=file, width = 4.5, height = 5, pointsize = 12, bg = "white")
opar <- par(mar=c(3,3,1,1)+0.1)
@

<<label=zfig,echo=FALSE,eval=FALSE>>= 
par(opar)
dev.null <- dev.off()
cat("\\includegraphics{", file, "}\n\n", sep="")
@

\section{Introduction}

Data on child mortalities for 167 census districts have been collected for a nine year period. Some census districts had very few cases, others larger numbers of cases. The underlying population at risk is taken as the 1981 census totals of children under 5. These data are only recorded for a single year, so need to be altered by a factor of nine to make them comparable to the mortality data. The data set is used by Marshall (1991), and presented in detail by Bailey and Gatrell (1995). It is their analysis we will be following here.

\begin{footnotesize}
<<echo=TRUE>>= 
library(maptools)
library(spdep) 
data(auckland)
@ 
\end{footnotesize}

\section{Analysis of counts and proportions}

If the counts data were collected from census aggregates all having the same population of under-five children, the mortality rates (numbers per thousand per year) would already be standardised for variations in the sizes of the populations at risk. Since they are not --- as we can see from the summary of \code{Under.5.1981} in the \code{auckland} data frame, we need to standardise.


\begin{footnotesize}
<<echo=TRUE>>= 
summary(auckland$Under.5.1981)
summary(auckland$Deaths.1977.85)
res <- probmap(auckland$Deaths.1977.85, 9*auckland$Under.5.1981)
str(res)
@ 
\end{footnotesize}
The \code{probmap()} function provides four measures as columns of a data frame, when given a vector of counts $y_i$, and a vector of population at risk levels $n_i$. The first that we map (Figure \ref{raw}), using ``pretty'' class intervals, is for crude (raw) rates $r_i = y_i / n_i$, per thousand per year. 

\begin{footnotesize}
<<echo=TRUE>>= 
brks <- c(-Inf, 2, 2.5, 3, 3.5, Inf)
cols <- c("#F7F7F7", "#CCCCCC", "#969696", "#636363", "#252525")
@ 
<<label=raw,fig=FALSE,echo=TRUE,eval=FALSE>>= 
plot(auckpolys, col=cols[findInterval(res$raw*1000, brks)], forcefill=FALSE)
legend(c(70,90), c(70,95), fill=cols, legend=leglabs(brks), bty="n")
@ 
\end{footnotesize}

\begin{figure}[htbp]
\begin{center} 
<<results=tex,echo=FALSE>>= 
<<afig>>
<<raw>>
<<zfig>>
@ 
\end{center}
\caption{Crude (raw) estimates of infant mortality per 1000 per year, 1977--1985, Auckland}
\label{raw}
\end{figure}


\begin{footnotesize}
<<echo=TRUE>>= 
brks <- c(-Inf,47,83,118,154,190,Inf)
cols <- c("#F7F7F7", "#CCCCCC", "#969696", "#636363", "#252525")
@ 
<<label=rel,fig=FALSE,echo=TRUE,eval=FALSE>>= 
plot(auckpolys, col=cols[findInterval(res$relRisk, brks)], forcefill=FALSE)
legend(c(70,90), c(70,95), fill=cols, legend=leglabs(brks), bty="n")
@ 
\end{footnotesize}
Next, we calculate the mean (expected) count as $\hat{\mu}_i = n_i ( \frac{\sum{y_i}}{\sum{n_i}})$, (Bailey and Gatrell p. 300). These values are also returned by the function, but not mapped here. Dividing the observed count in each census district by the expected count, and multiplying by 100, we get the standardised mortality rate, shown in Figure \ref{rel}; here the class intervals follow the book, Figure 8.1.

\begin{figure}[htbp]
\begin{center} 
<<results=tex,echo=FALSE>>= 
<<afig>>
<<rel>>
<<zfig>>
@ 
\end{center}
\caption{Standardised mortality ratios for child deaths per year, 1977--1985, Auckland}
\label{rel}
\end{figure}


\begin{footnotesize}
<<echo=TRUE>>= 
brks <- c(0,0.05,0.1,0.2,0.8,0.9,0.95,1)
cols <- c("#2166AC", "#67A9CF", "#D1E5F0", "#F7F7F7", "#FDDBC7", "#EF8A62",
 "#B2182B")
@ 
<<label=pmap,fig=FALSE,echo=TRUE,eval=FALSE>>= 
plot(auckpolys, col=cols[findInterval(res$pmap, brks)], forcefill=FALSE)
legend(c(70,90), c(70,95), fill=cols, legend=leglabs(brks), bty="n")
@ 
\end{footnotesize}
Finally, we show a probability map in Figure \ref{pmap}. This is derived by passing the observed and expected counts as arguments to the \code{ppois()} distribution function for the Poisson distribution, and has two tails --- this differs from the derivation given by Bailey and Gatrell on p. 302. The values can be read as indicating how probable it is that a value as large as the observed count compared to the expected count could have arisen at random. Values close to zero indicate unusually small observed counts in relation to the expected counts, values close to unity indicate unusually large observed counts in relation to the expected counts.

\begin{figure}[htbp]
\begin{center} 
<<results=tex,echo=FALSE>>= 
<<afig>>
<<pmap>>
<<zfig>>
@ 
\end{center}
\caption{Poisson probabilities for child deaths per year, 1977--1985, Auckland}
\label{pmap}
\end{figure}

\section{Empirical Bayes estimates}

Section 8.2.2, pp. 303--308, describes in detail the background for calculating 
Empirical Bayes estimates of rates. Briefly, in census districts with small populations at risk, we place less weight on the local rate, and more on the global rate for all the districts considered together. 

\begin{footnotesize}
<<echo=TRUE>>= 
res <- EBest(auckland$Deaths.1977.85, 9*auckland$Under.5.1981)
@ 
\end{footnotesize}
Using the method of moments, the estimated values of the mean $\hat{\phi} = $ \Sexpr{formatC(attr(res, "parameters")$a, 4, 3, "e")}, and variance $\hat{\gamma} = $ \Sexpr{formatC(attr(res, "parameters")$b, 4, 3, "e")} (returned as attributes of the data frame yielded by \code{EBest()}) correspond to those on p. 306. The mapped values (Figure \ref{EBest}) are essentially the same as Figure 8.3, but with the grey scale reversed, ``shrinking'' the rates of districts with small populations at risk towards the overall mean. 


\begin{footnotesize}
<<echo=TRUE>>= 
brks <- c(-Inf,2,2.5,3,3.5,Inf)
cols <- c("#F7F7F7", "#CCCCCC", "#969696", "#636363", "#252525")
@ 
<<label=EBest,fig=FALSE,echo=TRUE,eval=FALSE>>= 
plot(auckpolys, col=cols[findInterval(res$estmm*1000, brks)], forcefill=FALSE)
legend(c(70,90), c(70,95), fill=cols, legend=leglabs(brks), bty="n")
@
\end{footnotesize}

\begin{figure}[htbp]
\begin{center} 
<<results=tex,echo=FALSE>>= 
<<afig>>
<<EBest>>
<<zfig>>
@ 
\end{center}
\caption{Empirical Bayes global moment estimator of infant mortality per year, 1977--1985, Auckland}
\label{EBest}
\end{figure}


\section{Local Empirical Bayes estimates}

It is also possible, as Bailey and Gatrell indicate on p. 307, and in the exercises to the chapter on pp. 328--330, to ``shrink'' the estimates of census district rates not to the overall mean, but to neighbourhood means, where neighbourhoods are taken to be defined by districts sharing a boundary with the district of interest. Here, these shared boundaries have been listed in advance in the neighbours list object \code{auckland.nb}.

\begin{footnotesize}
<<echo=TRUE>>= 
summary(auckland.nb)
res <- EBlocal(spNamedVec("Deaths.1977.85", auckland),
 9*spNamedVec("Under.5.1981", auckland), auckland.nb)
@
\end{footnotesize}
As we can see, some neighbourhoods contain few districts, others relatively many. The idea is to ``shrink'' towards mean values that are appropriate for the part of the city that the district lies in, rather than to the overall mean value. On the other hand, the results will critically depend on our definition of neighbourhood.


\begin{footnotesize}
<<echo=TRUE>>= 
brks <- c(-Inf,2,2.5,3,3.5,Inf)
cols <- c("#F7F7F7", "#CCCCCC", "#969696", "#636363", "#252525")
@ 
<<label=EBloc,fig=FALSE,echo=TRUE,eval=FALSE>>= 
plot(auckpolys, col=cols[findInterval(res$est*1000, brks)], forcefill=FALSE)
legend(c(70,90), c(70,95), fill=cols, legend=leglabs(brks), bty="n")
@
\end{footnotesize}
Figure \ref{EBloc} shows the results of the use of this function for the given definition of neighbours. You may like to compare the two Empirical Bayes estimates of the rates with the crude rates shown in Figure \ref{raw}. Do you still feel as confident as you did in the beginning about which census districts are reporting high rates?

\begin{figure}[htbp]
\begin{center} 
<<results=tex,echo=FALSE>>= 
<<afig>>
<<EBloc>>
<<zfig>>
@ 
\end{center}
\caption{Empirical Bayes local moment estimator of infant mortality per year, 1977--1985, Auckland}
\label{EBloc}
\end{figure}

\section*{References}

\begin{description}

\item Marshall R M (1991) Mapping disease and mortality rates using
Empirical Bayes Estimators, \emph{Applied Statistics}, 40, 283--294.

\item Bailey T, Gatrell A (1995) \emph{Interactive Spatial Data Analysis},
Harlow: Longman.

\end{description}


<<echo=FALSE>>= 
options("width"=owidth)
options("warn"=ow)
@
\end{document}


