\documentclass[11pt]{article}

% \VignetteIndexEntry{Analysing Replicated Point Patterns in Spatstat}

\usepackage{graphicx}
\usepackage{Sweave}
\usepackage{bm}
\usepackage{anysize}

\marginsize{2cm}{2cm}{2cm}{2cm}

\newcommand{\pkg}[1]{\texttt{#1}}
\newcommand{\code}[1]{\texttt{#1}}
\newcommand{\R}{{\sf R}}
\newcommand{\spst}{\pkg{spatstat}}
\newcommand{\Spst}{\pkg{Spatstat}}

\newcommand{\bold}[1]{{\textbf {#1}}}

\newcommand{\indicate}[1]{\boldmaths{1}\{ {#1} \}}
\newcommand{\dee}[1]{\, {\rm d}{#1}}
\newcommand{\boldmaths}[1]{{\ensuremath\boldsymbol{#1}}}
\newcommand{\xx}{\boldmaths{x}}

\begin{document}
\bibliographystyle{plain}
\thispagestyle{empty}

<<echo=FALSE,results=hide,fig=FALSE>>=
options(SweaveHooks=list(fig=function() par(mar=c(1,1,1,1))))
@ 
\SweaveOpts{eps=TRUE}
\setkeys{Gin}{width=0.6\textwidth}

<<echo=FALSE,results=hide>>=
library(spatstat)
spatstat.options(image.colfun=function(n) { grey(seq(0,1,length=n)) })
sdate <- read.dcf(file = system.file("DESCRIPTION", package = "spatstat"),
         fields = "Date")
sversion <- read.dcf(file = system.file("DESCRIPTION", package = "spatstat"),
         fields = "Version")
options(useFancyQuotes=FALSE)
@ 

\title{Analysing replicated point patterns in \texttt{spatstat}}
\author{Adrian Baddeley}
\date{For \spst\ version \texttt{\Sexpr{sversion}}}
\maketitle

\begin{abstract}
This document describes \spst's capabilities for
fitting models to replicated point patterns.
More generally it applies to data from a designed experiment
in which the response from each unit is a spatial point pattern.
\end{abstract}

\tableofcontents
\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

`Replicated point patterns' are datasets consisting of 
several point patterns which can be 
regarded as independent repetitions of the same experiment. For example,
three point patterns taken from micrographs of three pipette samples of the 
same jug of milk, could be assumed to be replicated observations.

More generally we could have several experimental groups, with 
replicated point pattern data in each group. For example there may be
two jugs of milk that were treated differently, and we take three
pipette samples from each jug.

Even more generally our point patterns could be the result
of a designed experiment involving
control and treatment groups, covariates such as temperature,
and even spatial covariates (such as image data). 

This document describes some capabilities available in the \spst\ package
for analysing such data. 
\textbf{For further detail, see Chapter 16 of the spatstat book \cite{TheBook}.}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Overview of software}

The main components needed are:

\begin{itemize}
\item the model-fitting function \texttt{mppm}, an extension of the
  \texttt{spatstat} function \texttt{ppm}, that will fit Gibbs point process
  models to multiple point pattern datasets;
\item support for the class \texttt{"mppm"} of point process models 
   fitted by \texttt{mppm} (e.g. functions to print and plot the fitted model,
   analysis of deviance for Poisson models)
\item some tools for exploratory data analysis;
\item basic support for the data from such experiments
  by storing the data in a \emph{``hyperframe''}. A hyperframe is like
  a data frame, except that each entry in a column can be a point pattern
  or a pixel image, as well as a single number or categorical value.
\item four example datasets.
\end{itemize}

\section{Formulating the problem}

We view the experiment as involving a series of
{\em `units'\/}.
Each unit is subjected to a known set of experimental conditions 
(described by the values of the {\em covariates\/}), and
each unit yields a {\em response\/} which is a spatial point pattern.
The value of a particular covariate for each unit can be
either a single value (numerical, logical or factor),
or a pixel image.

Three important cases are:
\begin{description}
\item[independent replicates:]
We observe $n$
different point patterns that can be regarded as independent replicates,
i.e.\ independent realisations of the same point process.
The `responses' are the point patterns; there are no covariates.
\item[replication in groups:]
there are $K$ different experimental groups (e.g. control, aspirin,
nurofen). In group $k$ ($k=1,\ldots,K$) we observe $n_k$
point patterns which can be regarded as independent replicates within
this group. We regard this as an experiment with $n = \sum_k n_k$ 
units. The responses are the point patterns; there is one covariate
which is a factor (categorical variable) identifying which group 
each point pattern belongs to.
\item[general case:] 
there are covariates other than factors that influence
the response. The point patterns are assumed to be independent, 
but no two patterns have the same distribution.
\end{description}

Examples of these three cases are given in the 
datasets \texttt{waterstriders}, \texttt{pyramidal} and \texttt{demohyper}
respectively, which are installed in \spst.

\section{Installed datasets}

The following datasets are currently installed in \spst.

\begin{itemize}
\item \texttt{waterstriders}: Penttinen's \cite{pent84} waterstriders data
recording the locations of insect larvae on a pond in 3 independent
experiments.
\item \texttt{pyramidal}: data from Diggle, Lange and Benes 
\cite{digglangbene91} on the locations of pyramidal neurons in 
human brain, 31 human subjects grouped into 3 groups (controls,
schizoaffective and schizophrenic).
\item \texttt{flu}: data from Chen et al \cite{chenetal08}
giving the locations of two different virus proteins 
on the membranes of cells infected with influenza virus;
41 multitype point patterns divided into two virus types
(wild and mutant) and two stain types. 
\item \texttt{simba}: simulated data from an experiment with two groups
and 5 replicate point patterns per group.
\item \texttt{demohyper}: simulated data from an experiment with two 
  groups in which each experimental unit has a point pattern response
  and a pixel image covariate. 
\end{itemize}

\section{Lists of point patterns}

First we need a convenient way to store the \emph{responses}
from all the units in an experiment.

An individual point pattern is stored as an object of class \verb!"ppp"!. 
The easiest way to store all the responses is to form a list
of \verb!"ppp"! objects.

\subsection{Waterstriders data}

The \texttt{waterstriders} data are an example of this type.
The data consist of 3 independent point patterns representing the 
locations of insect larvae on a pond. See \texttt{help(waterstriders)}.

<<>>=
waterstriders
@ 

The \texttt{waterstriders} dataset is a list of point patterns.
It is a list, each of whose entries is a point pattern (object of class
\verb!"ppp"!). Note that the observation windows of the
three point patterns are {\tt not\/} identical. 

\subsection{The class \texttt{listof}}

For convenience, the \texttt{waterstriders} dataset also belongs to the
class \verb!"listof"!. This is a simple mechanism to allow us to
handle the list neatly --- for example, we can provide
special methods for printing, plotting and summarising the list.

\SweaveOpts{width=6,height=2}
\setkeys{Gin}{width=0.9\textwidth}

<<fig=TRUE>>=
plot(waterstriders, main="")
@ 

Notice that the 
plot method displays each entry of the list in a separate panel.
There's also the summary method:

<<>>=
summary(waterstriders)
@ 

\subsection{Creating a \texttt{listof} object}

For example, here is a simulated dataset containing three
independent realisations of the Poisson process with intensity 100.

<<>>=
X <- listof(rpoispp(100), rpoispp(100), rpoispp(100))
@ 

Then it can be printed and plotted.

<<fig=TRUE>>=
plot(X)
X
@ 

To convert an existing list to the class \code{listof}, use
\code{as.listof}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Hyperframes}

A \emph{hyperframe} is like a data frame, except that its entries
can be objects of any kind.
A hyperframe is effectively a two-dimensional array
in which each column consists of
values of one type (as in a data frame) or consists of
objects of one class. 

The entries in a hyperframe can be point patterns, pixel images,
windows, or any other objects. 

To analyse an experiment, we will store {\bf all} the data from the experiment
in a single hyperframe. The rows of the hyperframe will correspond 
to different experimental units,
while the columns represent different variables 
(response variables or covariates).

\subsection{Creating hyperframes}

The function \texttt{hyperframe} will create a hyperframe.

<<eval=FALSE>>=
hyperframe(...)
@ 

The arguments \verb!...! are any number of arguments of
the form \texttt{tag=value}. Each \texttt{value} will
become a column of the array. The \texttt{tag} determines the name
of the column.

Each \texttt{value} can be either
\begin{itemize}
\item an atomic vector or factor
  (i.e. numeric vector, integer vector, character vector, logical
  vector, complex vector or factor)
\item a list of objects which are all of the same class
\item one atomic value, which will be replicated to make an atomic
  vector or factor
\item one object, which will be replicated to make a list of identical objects.
\end{itemize}
    
All columns (vectors, factors and lists) must be of the same length,
if their length is greater than 1. 

For example, here is a hyperframe containing a column of
numbers and a column of \emph{functions}:

<<>>=
H <- hyperframe(X=1:3, Y=list(sin,cos,tan))
H
@ 

Note that a column of character strings will be converted to a factor, 
unless you set \texttt{stringsAsFactors=FALSE} in the call to
\code{hyperframe}. This is the same behaviour as for the function 
\code{data.frame}.

<<>>=
G <- hyperframe(X=1:3, Y=letters[1:3], Z=factor(letters[1:3]),
                W=list(rpoispp(100),rpoispp(100), rpoispp(100)),
                U=42,
                V=rpoispp(100), stringsAsFactors=FALSE)
G
@ 

This hyperframe has 3 rows and 6 columns. 
The columns named \texttt{U} and \texttt{V}
are constant (all entries in a column are the same). The column named
\texttt{Y} is a character vector while \texttt{Z} is a factor.

\subsection{Hyperframes of data}

To analyse an experiment, we will store {\bf all} the data from the experiment
in a single hyperframe. The rows of the hyperframe will correspond 
to different experimental units,
while the columns represent different variables 
(response variables or covariates).

Several examples of hyperframes are provided with the package,
including \texttt{demohyper}, \texttt{flu}, \texttt{simba}
and \texttt{pyramidal}, described above.

The \texttt{simba} dataset contains simulated data from an
experiment with a `control' group and a `treatment' group, each 
group containing 5 experimental units. The responses in the control group are 
independent Poisson point patterns with intensity 80.
The responses in the treatment group are independent realisations of
a Strauss process (see \texttt{help(simba)} for details).
The \texttt{simba} dataset is a hyperframe with 10 rows and 2 columns:
\texttt{Points} (the point patterns) and \texttt{group} (a factor 
with levels \texttt{control} and \texttt{treatment}).

<<>>=
simba
@ 

The \texttt{pyramidal} dataset contains data from Diggle, Lange and Benes 
\cite{digglangbene91} on the locations of pyramidal neurons in 
human brain. One point pattern was observed in each of 31 human subjects.
The subjects were classified
into 3 groups (controls, schizoaffective and schizophrenic).
The \texttt{pyramidal} dataset is a hyperframe with 31 rows
and 2 columns: \code{Neurons} (the point patterns) and \code{group}
(a factor with levels \texttt{control}, \texttt{schizoaffective} 
and \texttt{schizophrenic}).

<<>>=
pyramidal
@ 

The \texttt{waterstriders} dataset is not a hyperframe; it's just a 
list of point patterns. It can easily be converted into a hyperframe:

<<>>=
ws <- hyperframe(Striders=waterstriders)
@ 

\subsection{Columns of a hyperframe}

Individual columns of a hyperframe can be extracted using \verb!$!:

<<>>=
H$X
H$Y
@ 

The result of \verb!$! is a vector or factor if the column contains
atomic values; otherwise it is a list of objects (with class \texttt{"listof"}
to make it easier to print and plot).

Individual columns can also be assigned (overwritten or created)
using \verb!$<-!:

<<>>=
H$U <- letters[1:3]
H
@ 

This can be used to build up a hyperframe column-by-column:

<<>>=
G <- hyperframe()
G$X <- waterstriders
G$Y <- 1:3
G
@ 

\subsection{Subsets of a hyperframe}

Other subsets of a hyperframe
can be extracted with \verb![!:

<<>>=
H[,1]
H[2,]
H[2:3, ]
H[1,1]
@ 

The result of \verb![! is a hyperframe, unless you set \verb!drop=TRUE!
and the subset consists of only one element or one column:

<<>>=
H[,1,drop=TRUE]
H[1,1,drop=TRUE]
H[1,2,drop=TRUE]
@ 

There is also a method for \verb![<-! that allows
you to assign values to a subset of a hyperframe.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Plotting}

\subsection{Plotting a \code{listof} object}

The plot method for \code{listof} objects has formal arguments

<<eval=FALSE>>=
plot.listof(x, ..., main, arrange = TRUE, nrows = NULL, ncols = NULL)
@ 

where \code{main} is a title for the entire page.

If \code{arrange=TRUE} then the entries of the list are displayed  
in separate panels on the same page (with \code{nrows} rows and
\code{ncols} columns of panels), while if \code{arrange=FALSE} then 
the entries are just plotted as a series of plot frames.

The extra arguments \verb!...! control the individual plot panels.
These arguments will be passed to the plot method
that displays each entry of the list. Suitable arguments depend on the
type of entries.

<<fig=TRUE>>= 
plot(waterstriders, pch=16, nrows=1)
@ 

\subsection{Plotting a hyperframe}

\subsubsection{Plotting one column}

If \code{h} is a hyperframe, then the default action of
\code{plot(h)} is to extract the first column
of \code{h} and plot each of the entries in a separate panel on 
one page (actually using the plot method for class \verb!"listof"!). 

\SweaveOpts{width=7,height=5}
\setkeys{Gin}{width=0.9\textwidth}

<<fig=TRUE>>=
plot(simba)
@ 

This only works if the entries in the first column are objects 
for which a plot method is defined (for example, point patterns, images,
windows). 

To select a different column, use \verb!$! or \verb![!:

\SweaveOpts{width=6,height=2}
\setkeys{Gin}{width=0.9\textwidth}

<<fig=TRUE>>=
H <- hyperframe(X=1:3, Y=list(sin,cos,tan))
plot(H$Y)
@ 

The plot can be controlled using the arguments for \code{plot.listof}
(and, in this case, \code{plot.function}, since \verb!H$Y! consists of 
functions).

\subsubsection{Complex plots}

More generally, we can display any kind of higher-order plot
involving one or more columns of a hyperframe:

<<eval=FALSE>>=
plot(h, e)
@ 

where \code{h} is a hyperframe and \code{e} is an \R\ language call 
or expression that must be evaluated in each row to generate each plot panel.

\SweaveOpts{width=9,height=5}
\setkeys{Gin}{width=0.9\textwidth}

<<fig=TRUE>>=
plot(demohyper, quote({ plot(Image, main=""); plot(Points, add=TRUE) }))
@ 

Note the use of \code{quote}, which prevents the code
inside the braces from being evaluated immediately.

To plot the $K$-functions of each of the patterns in the
\code{waterstriders} dataset,

\SweaveOpts{width=6,height=2}
\setkeys{Gin}{width=0.9\textwidth}

<<fig=TRUE>>=
H <- hyperframe(Bugs=waterstriders)
plot(H, quote(plot(Kest(Bugs))), marsize=1)
@ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Data analysis}

\subsection{Computing with hyperframes}

Often we want to perform some computation on each row of
a hyperframe. 

In a data frame, this can be done using the command \code{with}:

<<>>=
df <- data.frame(A=1:10, B=10:1)
with(df, A-B)
@ 

In this example, the expression \code{A-B} is evaluated
in each row of the data frame, and the result is a vector
containing the computed values for each row.
The function \code{with} is generic, and has a method for data frames,
\code{with.data.frame}. The computation above
was executed by \code{with.data.frame}. 

The same syntax is available for hyperframes
using the method \code{with.hyperframe}:

<<eval=FALSE>>=
with(h,e)
@ 

Here \code{h} is a hyperframe, and 
\code{e} is an {\sf R} language construct involving the names
of columns in \code{h}.
For each row of \code{h}, the expression \code{e} will be evaluated
in such a way that each entry in the row is identified by its
column name. 

<<>>=
H <- hyperframe(Bugs=waterstriders)
with(H, npoints(Bugs))
with(H, distmap(Bugs))
@

The result of \code{with.hyperframe} 
is a list of objects (of class \verb!"listof"!),
or a vector or factor if appropriate.

Notice that (unlike the situation for data frames)
the operations in the expression \code{e} do not have to
be vectorised. For example, \code{distmap} expects 
a single point pattern, and is not vectorised to deal with 
a list of point patterns. Instead, the expression \code{distmap(Bugs)}
is evaluated separately in each row of the hyperframe.

\subsection{Summary statistics}

One application of \code{with.hyperframe} is to calculate summary statistics
for each row of a hyperframe.

For example, the number of points in a point pattern \code{X}
is returned by \code{npoints(X)}. To calculate this for each of the 
responses in the \code{simba} dataset,

<<>>=
with(simba, npoints(Points))
@ 

The summary statistic can be any kind of object. For example, to
compute the empirical $K$-functions for each of the 
patterns in the \code{waterstriders} dataset,

<<>>=
H <- hyperframe(Bugs=waterstriders)
K <- with(H, Kest(Bugs))
@ 

To plot these $K$-functions you can then just type 

\SweaveOpts{width=6,height=2}
\setkeys{Gin}{width=0.9\textwidth}

<<fig=TRUE>>=
plot(K)
@ 

The summary statistic for each row could be a numeric vector:

<<>>=
H <- hyperframe(Bugs=waterstriders)
with(H, nndist(Bugs))
@ 

The result is a list, each entry being a vector of nearest neighbour distances.
To find the minimum interpoint distance in each pattern:

<<>>=
with(H, min(nndist(Bugs)))
@ 

\subsection{Generating new columns}

New columns of a hyperframe can be created by computation
from the existing columns. 

For example, I can add a new column to the \code{simba} dataset
that contains pixel images of the distance maps for each of the
point pattern responses. 

<<fig=FALSE>>=
simba$Dist <- with(simba, distmap(Points))
@ 

\subsection{Simulation}


This can be useful for simulation. For example, to generate
Poisson point patterns with different intensities, where the 
intensities are given by a numeric vector \code{lambda}:

\SweaveOpts{width=6,height=6}
\setkeys{Gin}{width=0.7\textwidth}

<<fig=TRUE>>=
lambda <- rexp(6, rate=1/50)
H <- hyperframe(lambda=lambda)
H$Points <- with(H, rpoispp(lambda))
plot(H, quote(plot(Points, main=paste("lambda=", signif(lambda, 4)))))
@ 

It's even simpler to generate 10 independent Poisson point patterns
with the \emph{same} intensity 50, say:

<<fig=FALSE>>=
H$X <- with(H, rpoispp(50))
@ 

\noindent The expression \code{rpoispp(50)} is evaluated once in each row,
yielding a different point pattern in each row because of the
randomness.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Exploratory data analysis}

Before fitting models to the data, it is prudent to explore 
the data to detect unusual features and to suggest appropriate 
models. 

\subsection{Exploring spatial trend and covariate effects}

Points may be distributed non-uniformly either because they are
intrinsically non-uniform (``spatial trend'') or because their abundance
depends on a spatial covariate (``covariate effects''). 

Non-uniformity of a point pattern can be investigated using 
the kernel smoothed intensity. This is the convolution of the point pattern
with a smooth density called the kernel. Effectively each point
in the pattern is replaced by a copy of the kernel, and the sum of all
copies of the kernel is the kernel-smoothed intensity function.
It is computed by \texttt{density.ppp} separately for each point pattern.

<<fig=TRUE>>=
plot(simba, quote(plot(density(Points), main="")), nrows=2)
@ 

Covariate effects due to a real-valued spatial covariate (a real-valued
pixel image) can be investigated
using the command \code{rhohat}. This uses a kernel smoothing
technique to fit a model of the form 
\[
           \lambda(u) = \rho(Z(u))
\]
where $\lambda(u)$ is the point process intensity at a location $u$,
and $Z(u)$ is the value of the spatial covariate at that location.
Here $\rho$ is an unknown, smooth function which is to be estimated.
The function $\rho$ expresses the effect of the
spatial covariate on the point process intensity. If $\rho$ turns out to
be constant, then the covariate has no effect on point process intensity
(and the constant value of $\rho$ is the constant intensity of the
point process). 

<<fig=TRUE>>=
rhos <- with(demohyper, rhohat(Points, Image))
plot(rhos)
@ 

\SweaveOpts{width=6,height=4}
\setkeys{Gin}{width=0.9\textwidth}

\subsection{Exploring interpoint interaction}

Still to be written.
See Chapter 16 of \cite{baddrubaturn15}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Fitting models of spatial trend}

The command \code{mppm} fits models to multiple point patterns.
Its syntax is very similar to that of \code{lm} and \code{glm}:

<<eval=FALSE>>=
mppm(formula, data, interaction, ...)
@ 

where \code{formula} is a formula describing the systematic trend
part of the model, \code{data} is a hyperframe containing all the data 
(responses and covariates), and \code{interaction} determines the
stochastic interpoint interaction part of the model.

For example:

<<eval=FALSE>>=
mppm(Points ~ group, simba, Poisson())
@ 

Note that the formula has a left hand side, which identifies
the response. This should be the name of a column of \code{data}.

\subsection{Trend formula}

The right side of \code{formula} is an expression for the
linear predictor (effectively the {\bf logarithm} 
of the spatial trend). 

The variables appearing in the right hand side
of \code{formula} should be either
\begin{itemize}
\item names of columns in \code{data}
\item objects in the {\sf R} global
environment (such as \code{pi} and \code{log})
\item the reserved names \code{x}, \code{y} 
(representing Cartesian coordinates), \code{marks} (representing mark values
attached to points) or \code{id} (a factor representing the row number
in the hyperframe).
\end{itemize}

\subsubsection{Design covariates}

The variables in the trend could be `design covariates'.

For example, to fit a model to the \code{simba} dataset
in which all patterns are independent replicates of the
same uniform Poisson process, with the same constant
intensity:

<<>>=
mppm(Points ~ 1, simba)
@ 

To fit a model in which the two groups of patterns (control and treatment
groups) each consist of independent replicates of a uniform Poisson process,
but with possibly different intensity in each group:

<<>>=
mppm(Points ~ group, simba)
@ 

To fit a uniform Poisson process to each pattern, with 
different intensity for each pattern:

<<>>=
mppm(Points ~ id, simba)
@ 

\subsubsection{Spatial covariates}

The variables in the trend could be `spatial covariates'.

For example, the \code{demohyper} dataset has a column \code{Image}
containing pixel images. 

<<>>=
mppm(Points ~ Image, data=demohyper)
@ 

This model postulates that each pattern is a Poisson process
with intensity of the form
\[
      \lambda(u) = \exp(\beta_0 + \beta_1 Z(u))
\]
at location $u$, where $\beta_0, \beta_1$ are coefficients
to be estimated, and $Z(u)$ is the value of the pixel image
\code{Image} at location $u$.

It may or may not be appropriate to assume that the intensity of the points
is an exponential function of the image pixel value $Z$.
If instead 
we wanted the intensity $\lambda(u)$ to be \emph{proportional} to $Z(u)$,
the appropriate model is 

<<eval=FALSE>>=
mppm(Points ~ offset(log(Image)), data=demohyper)
@ 

which corresponds to an intensity proportional to \code{Image},
\[
      \lambda(u) = \exp(\beta_0 + \log Z(u)) = e^{\beta_0} \; Z(u).
\]
The \code{offset} indicates that there is no coefficient in front
of $\log Z(u)$. 

Alternatively we could allow a coefficient:

<<eval=FALSE>>=
mppm(Points ~ log(Image), data=demop)
@ 

which corresponds to a gamma transformation of \code{Image},
\[
      \lambda(u) = \exp(\beta_0 + \beta_1 \log Z(u)) 
      = e^{\beta_0} \; Z(u)^{\beta_1}.
\]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Interpoint interaction}

The stochastic interpoint interaction in a point process model
is specified by the arguments \code{interaction} and (optionally)
\code{iformula} in 

<<eval=FALSE>>=
mppm(formula, data, interaction, ..., iformula=NULL)
@ 

\subsection{Same interaction for all patterns}

In the simplest case,
the argument \texttt{interaction} is one of the familiar 
objects that describe the point process interaction structure.
It is an object of class \texttt{"interact"}
created by calling one of the functions

\begin{center}
  \begin{tabular}{rl}
    \texttt{Poisson()}  & the Poisson point process\\     
    \texttt{Hardcore()}  & the hard core process \\ 
    \texttt{Strauss()}  & the Strauss process \\ 
    \texttt{StraussHard()}  & the Strauss/hard core point process\\ 
    \texttt{Softcore()}  & pairwise interaction, soft core potential\\  
    \texttt{PairPiece()}  & pairwise interaction, piecewise constant \\ 
    \texttt{DiggleGatesStibbard() }  & Diggle-Gates-Stibbard pair potential \\ 
    \texttt{DiggleGratton() }  & Diggle-Gratton pair potential \\ 
    \texttt{Fiksel() }  & Fiksel pair potential \\ 
    \texttt{LennardJones() }  & Lennard-Jones pair potential \\ 
    \texttt{Pairwise()}  &	pairwise interaction, user-supplied potential\\ 
    \texttt{AreaInter()}  &		area-interaction potential\\ 
    \texttt{Geyer()}	  & Geyer's saturation process\\ 
    \texttt{BadGey()}	  & multiscale Geyer saturation process\\ 
    \texttt{Saturated()}  &	Saturated pair model, user-supplied potential\\ 
    \texttt{OrdThresh()}  &		Ord process, threshold potential\\ 
    \texttt{Ord()}  & 		        Ord model, user-supplied potential \\ 
    \texttt{MultiStrauss()}  & 		multitype Strauss process \\ 
    \texttt{MultiStraussHard()}  & 	multitype Strauss/hard core process \\ 
    \texttt{Concom()}  &	connected component interaction \\ 
    \texttt{Hybrid()}  &	hybrid of several interactions \\ 
  \end{tabular}
\end{center}

In this `simple' usage of \texttt{mppm}, the point process model
assumes that all point patterns have exactly the same interpoint interaction,
(with the same interaction parameters), and only differ in their spatial trend.

\subsection{Hyperframe of interactions}

More generally the argument \code{interaction} can be a hyperframe
containing objects of class \texttt{"interact"}. 

For example, we
might want to fit a Strauss process to each point pattern, but with
a different Strauss interaction radius for each pattern. 

<<echo=FALSE,results=hide>>=
radii <- with(simba, mean(nndist(Points)))
@ 

Then \code{radii} is a vector of numbers which we could use
as the values of the interaction radius for each case. 
First we need to make the interaction objects:

<<>>=
Rad <- hyperframe(R=radii)
Str <- with(Rad, Strauss(R))
@ 

Then we put them into a hyperframe and fit the model:

<<>>=
Int <- hyperframe(str=Str)
mppm(Points ~ 1, simba, interaction=Int)
@ 

An important constraint is that all of the interaction objects 
in one column must be \emph{instances of the same process} (e.g. Strauss)
albeit possibly having different parameter values. For example, you cannot
put Poisson and Strauss processes in the same column. 

\subsection{Interaction formula}

If \code{interaction} is a hyperframe, then 
the additional argument \code{iformula} may be used to
fully specify the interaction.

(An \code{iformula} is also required if \code{interaction}
has more than one column.)

The \code{iformula} should be a formula without a left hand side.
Variables on the right hand side are typically the names of
columns in \code{interaction}.

\subsubsection{Selecting one column}

If the right hand side of \code{iformula} is a single name,
then this identifies the column in \code{interaction}
to be used as the interpoint interaction structure.

<<>>=
h <- hyperframe(Y=waterstriders)
g <- hyperframe(po=Poisson(), str4 = Strauss(4), str7= Strauss(7))
mppm(Y ~ 1, data=h, interaction=g, iformula=~str4)
@ 

\subsubsection{Interaction depending on design} 

The \code{iformula} can also involve columns of \code{data}, but
only those columns that are vectors or factors. This allows us to
specify an interaction that depends on the experimental design.
[This feature is {\bf experimental}.]
For example

<<>>=
fit <- mppm(Points ~ 1, simba, Strauss(0.07), iformula = ~Interaction*group)
@ 

Since \code{Strauss(0.1)} is not a hyperframe, it is first converted
to a hyperframe with a single column named \code{Interaction}.

The \code{iformula = ~Interaction*group} specifies (since \code{group}
is a factor) that the interpoint interaction shall have a different 
coefficient in each experimental group. That is, we fit a model
which has two different values for the Strauss interaction parameter $\gamma$,
one for the control group and one for the treatment group.

When you print the result of such a fit, 
the package tries to 
do `automatic interpretation' of the fitted model (translating the
fitted interaction coefficients into meaningful numbers like $\gamma$). 
This will be successful in \emph{most} cases:

<<>>=
fit
@ 

<<echo=FALSE,results=hide>>=
co <- coef(fit)
si <- function(x) { signif(x, 4) }
@ 

Thus we see that the estimate of the Strauss parameter $\gamma$ 
for the control group is \Sexpr{si(exp(co[2]))}, and 
for the treatment group \Sexpr{si(exp(sum(co[c(2,4)])))}
(the correct values in this simulated dataset were $1$ and $0.5$).

The fitted model can also be interpreted directly from the fitted
canonical coefficients:

<<>>=
coef(fit)
@ 

The last output shows all the coefficients $\beta_j$
in the linear predictor for the (log) conditional intensity. 

The interpretation of the model coefficients, for any fitted model in \R,
depends on the \emph{contrasts} which were applicable when the model was
fitted. This is part of the core {\sf R} system: see \code{help(contrasts)}
or \code{options(contrasts)}. If you did not specify otherwise, 
the default is to use \emph{treatment contrasts}. This means that,
for an explanatory variable which is a \texttt{factor} with $N$ levels, 
the first level of the factor is used as a baseline, and the
fitted model coefficients represent the factor levels $2, 3, \ldots, N$
relative to this baseline. 

In the output above, there is a coefficient for \code{(Intercept)} 
and one for \code{grouptreatment}. These are coefficients related to
the \code{group} factor. According to the ``treatment contrasts'' rule, 
the \code{(Intercept)} coefficient is
the estimated effect for the control group, and the
\code{grouptreatment} coefficient is the estimated difference between
the treatment and control groups. Thus the fitted first order
trend is $\exp(\Sexpr{si(co[1])}) = \Sexpr{si(exp(co[1]))}$ 
for the control group
and $\exp(\Sexpr{si(co[1])} + \Sexpr{si(co[3])}) 
 = \Sexpr{si(exp(sum(co[c(1,3)])))}$ for the treatment group.
The correct values in this simulated dataset were 
$80$ and $100$.

The remaining coefficients in the output are \code{Interaction}
and \code{Interaction:grouptreatment}. Recall that the Strauss process
interaction term 
is $\gamma^{t(u,\xx)} = \exp(t(u,\xx) \log\gamma)$
at a spatial location $u$, for a point pattern $\xx$. 
Since we're using treatment contrasts, the coefficient 
\code{Interaction} is the estimate of 
$\log\gamma$ for the control group.
The coefficient  \code{Interaction:grouptreatment} is the 
estimate of the difference in $\log\gamma$ between the 
treatment and control groups. Thus the estimated Strauss interaction
parameter $\gamma$ is $\exp(\Sexpr{si(co[2])}) = \Sexpr{si(exp(co[2]))}$ 
for the control group and 
$\exp(\Sexpr{si(co[2])} + (\Sexpr{si(co[4])})) = \Sexpr{si(exp(co[2]+co[4]))}$ 
for the treatment group.
The correct values were $1$ and $0.5$.

\subsubsection{Completely different interactions for different cases} 

In the previous example, when we fitted a Strauss model to all 
point patterns in the \code{simba} dataset, the fitted model for the
patterns in the control group was close to Poisson ($\gamma \approx 1$).
Suppose we now want to fit a model which {\it is}
Poisson in the control group, and Strauss in the treatment group. 
The Poisson and Strauss interactions must be given as separate columns
in a hyperframe of interactions:

<<eval=FALSE>>=
interaction=hyperframe(po=Poisson(), str=Strauss(0.07))
@ 

What do we write for the 
\code{iformula}? The following \emph{will not} work:
<<eval=FALSE>>=
iformula=~ifelse(group=="control", po, str)
@ 
This does not work because the Poisson and Strauss models are `incompatible'
inside such expressions. The canonical sufficient statistics for
the Poisson and Strauss processes do not have the same dimension. 
Internally in \code{mppm} we translate the symbols \code{po} and \code{str}
into matrices; the dimensions of these matrices are different,
so the \code{ifelse} expression cannot be evaluated.

Instead we need something like the following:
<<eval=FALSE>>=
iformula=~I((group=="control")*po) + I((group=="treatment") * str)
@ 
The letter \code{I} here is a standard R function that prevents its argument
from being interpreted as a formula (thus the \code{*} is interpreted
as multiplication instead of a model interaction). The expression
\code{(group=="control")} is logical, and when multiplied by the matrix
\code{po}, yields a matrix.

So the following does work:

<<>>=
g <- hyperframe(po=Poisson(), str=Strauss(0.07))
fit2 <- mppm(Points ~ 1, simba, g, 
             iformula=~I((group=="control")*po) 
                     + I((group=="treatment") * str))
fit2
@ 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%#%^!ifdef RANDOMEFFECTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Random effects}

\subsection{Mixed effects models}

It is also possible to fit models that include `random effects'.
Effectively, some of the coefficients
in the model are assumed to be Normally-distributed random variables
instead of constants.

\subsubsection{Mixed Poisson model}

Consider the simplest model of a uniform Poisson process which we fitted to
the 3 point patterns of waterstriders. It might be sensible to assume that
each pattern is a realisation of a Poisson process, but with 
{\em random intensity\/}. In each realisation the intensity $\lambda$
is constant across different locations, but it
is a different, random value in different realisations.
This example is called a `mixed Poisson process' and belongs to the
class of `Cox processes' (Poisson processes with random intensity 
functions).
Let's assume further that the log-intensity 
is a Normal random variable.
Then the model is a (very degenerate) special case of
a `log-Gaussian Cox process'. 

To fit such a model we use the standard techniques of mixed effects
models \cite{lairware82,davigilt95,pinhbate00}.
The mixed Poisson process which we discussed above would
be written in standard form
\begin{equation}
\label{mixPois}
   \lambda_i(u) = \exp(\mu + Z_i)
\end{equation}
for the $i$th point pattern, where $\mu$ is a parameter to be estimated 
(the `fixed effect')
and $Z_i \sim N(0, \sigma^2)$ is a zero-mean Normal random variable
(the `random effect' for point pattern $i$). In the simplest case we
would assume that $Z_1, \ldots, Z_n$ are independent.
The variance $\sigma^2$ of the random effects
would be estimated. One can also estimate the individual realised values
$z_i$ of the random effects for each point pattern, although these are
usually not of such great interest. Since the model includes both
fixed and random effects, it is called a ``mixed-effects'' model.

\subsubsection{Dependence structure}

When we formulate a random-effects or mixed-effects model, we must 
specify the dependence structure of the random effects. In the model above
we assumed that the $Z_i$ are independent for all point patterns $i$. 
If the experiment consists of two groups, we could alternatively assume
that $Z_i = Z_j$ whenever $i$ and $j$ belong to the same group. In other words
all the patterns in one group have the same value of the random effect.
So the random effect is associated with the group rather than with
individual patterns. This could be appropriate if, for example, 
the groups represent different
batches of a chemical. Each batch is prepared under slightly different 
conditions so we believe that there are random variations between batches,
but within a batch we believe that the chemical is well-mixed.

\subsubsection{Random effects are coefficients}

In the mixed Poisson model (\ref{mixPois}), 
the random effect is an additive constant
(with a random value) in the log-intensity. In general, a
random effect is a \emph{coefficient} of one of the covariates.
For example if $v$ is a real-valued design covariate (e.g. `temperature'), 
with value $v_i$ for the $i$th point pattern, then we could assume
\begin{equation}
\label{ranef2}
   \lambda_i(u) = \exp(\mu + Z_i v_i)
\end{equation}
where $Z_i \sim N(0, \sigma^2)$ are independent for different $i$.
This model has a random effect in the dependence on $v$.

We could also have a random effect for a spatial covariate $V$.
Suppose $V_i$ is a real-valued image for the $i$th pattern 
(so that $V_i(u)$ is the value of some covariate at the location $u$
for the $i$th case). Then we could assume
\begin{equation}
\label{ranef3}
   \lambda_i(u) = \exp(\mu + Z_i V_i(u))
\end{equation}
where $Z_i \sim N(0, \sigma^2)$ are independent for different $i$.
This kind of random effect would be appropriate if, for example, the 
images $V_i$ are not `normalised' or `standardised' relative to each 
other (e.g.\ they are images taken under different illumination). Then the
coefficients $Z_i$ effectively include the rescaling necessary to standardise
the images.

\subsection{Fitting a mixed-effects model}

The call to \texttt{mppm} can also include the argument
\texttt{random}. This should be a formula (with no left-hand side) 
describing the structure of random effects. 

The formula for random effects
must be recognisable to \texttt{lme}. It is typically of the form
\begin{verbatim}
      ~x1 + ... + xn | g
\end{verbatim}
or
\begin{verbatim}
      ~x1 + ... + xn | g1/.../gm
\end{verbatim}
where \verb!x1 + ... + xn! specifies the covariates for the random effects
and \texttt{g} or \verb!g1/.../gm! determines the grouping (dependence) 
structure. Here \code{g} or \code{g1, \ldots, gm} should be factors.

To fit the mixed Poisson model (\ref{mixPois}) to the waterstriders,
we want to have a random intercept coefficient
(so \texttt{x} is \texttt{1}) that varies for different point patterns
(so \texttt{g} is \texttt{id}). 
The reserved name \code{id} is a factor referring to the individual
point pattern. Thus

<<>>=
H <- hyperframe(P=waterstriders)
mppm(P ~ 1, H, random=~1|id)
@ 

To fit the mixed effects model (\ref{ranef2}) to the coculture data
with the \code{AstroIm} covariate, with a random effect associated
with each well,

<<eval=FALSE>>=
mppm(Neurons ~ AstroIm, random=~AstroIm|WellNumber)
@ 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%#%^!endif
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Studying the fitted model}

Fitted models produced by \code{mppm}
can be examined and validated in many ways.

\subsection{Fits for each pattern}


\subsubsection{Subfits}

The command \code{subfits} takes an \code{mppm} object
and extracts, for each individual point pattern, 
the fitted point process model for that pattern
\emph{that is implied by the overall fit}. It returns a list of 
objects of class \code{ppm}. 

<<>>=
H <- hyperframe(W=waterstriders)
fit <- mppm(W ~ 1, H)
subfits(fit)
@ 

In this example the result is a list of three \code{ppm} objects
representing the implied fits for each of the three point patterns
in the \code{waterstriders} dataset.
Notice that {\bf the fitted coefficients are the same} in all three 
models. 

Note that there are some unresolved difficulties with the implementation of
\code{subfits}. Two completely different implementations are supplied
in the package; they are called \code{subfits.old} 
%(used in versions 0.1--1 and earlier) 
and \code{subfits.new}.% (introduced in 0.1--2).
The old version would occasionally crash. 
Unfortunately the newer version \code{subfits.new} is quite memory-hungry
and sometimes causes R to hang.
We're still working on this problem. So for the time being,
\code{subfits} is the same as \code{subfits.old}. You can change this
simply by reassigning, e.g.

<<eval=FALSE>>=
subfits <- subfits.new
@ 

\subsubsection{Fitting separately to each pattern} 

For comparison, we could fit a point process model separately
to each point pattern dataset using \code{ppm}. The easy way to do this is 
with \code{with.hyperframe}.

To fit a \emph{separate} uniform Poisson 
point process to each of the three waterstriders patterns,

<<>>=
H <- hyperframe(W=waterstriders)
with(H, ppm(W))
@ 

The result is again a list of three fitted point process models
(objects of class \code{ppm}), but now the fitted coefficients
are different.

\subsection{Residuals}

One standard way to check a fitted model
is to examine  the residuals. 

\subsubsection{Point process residuals}

Some recent papers \cite{baddetal05,baddmollpake08} have defined
residuals for a fitted point process model (fitted to a \emph{single}
point pattern). These residuals are implemented in \code{spatstat}
as \code{residuals.ppm} and apply to an object of class \code{ppm}, 
that is, a model fitted to a \emph{single} point pattern. 

The command \code{residuals.mppm} computes the point process residuals
for an \code{mppm} object. 

<<>>=
fit <- mppm(P ~ x, hyperframe(P=waterstriders))
res <- residuals(fit)
@ 

The result is a list, with one entry for each of the point pattern
datasets. Each list entry contains the point process residuals 
for the corresponding point pattern dataset.
Each entry in the list is a signed measure (object of class \code{"msr"})
as explained in the help for \code{residuals.ppm}).
It can be plotted:

<<fig=TRUE>>=
plot(res)
@ 

You probably want the smoothed residual field:

<<fig=TRUE>>=
smor <- with(hyperframe(res=res), Smooth(res, sigma=4))
plot(smor)
@ 

\subsubsection{Sums of residuals}

It would be useful to have a residual that is a single value
for each point pattern (representing how much that point pattern
departs from the model fitted to all the point patterns).

That can be computed by \emph{integrating} the residual measures
using the function \code{integral.msr}:

<<>>=
fit <- mppm(P ~ x, hyperframe(P=waterstriders))
res <- residuals(fit)
totres <- sapply(res, integral.msr)
@ 

In designed experiments we can plot these total residuals against
the design covariates:

<<fig=TRUE>>=
fit <- mppm(Points~Image, data=demohyper)
resids <- residuals(fit, type="Pearson")
totres <- sapply(resids, integral.msr)
areas <- with(demohyper, area.owin(as.owin(Points)))
df <- as.data.frame(demohyper[, "Group"])
df$resids <- totres/areas
plot(resids~Group, df)
@ 

\subsubsection{Four-panel diagnostic plots} 

Sometimes a more useful tool is the function \code{diagnose.ppm}
which produces a four-panel diagnostic plot based on the 
point process residuals. However, it is only available for 
\code{ppm} objects.

To obtain a four-panel diagnostic plot for each of the 
point patterns, do the following:

\begin{enumerate}
\item fit a model to multiple point patterns using \code{mppm}.
\item extract the individual fits using \code{subfits}.
\item plot the residuals of the individual fits.
\end{enumerate}

For example:

<<fig=TRUE>>=
fit <- mppm(P ~ 1, hyperframe(P=waterstriders))
sub <- hyperframe(Model=subfits(fit))
plot(sub, quote(diagnose.ppm(Model)))
@ 

(One could also do this for models fitted separately to the 
individual point patterns.)

\subsubsection{Residuals of the parameter estimates}

We can also compare the parameter estimates obtained 
by fitting the model simultaneously to all patterns (using \code{mppm})
with those obtained by fitting the model separately to each 
pattern (using \code{ppm}). 

<<>>=
H <- hyperframe(P = waterstriders)
fitall <- mppm(P ~ 1, H)
together <- subfits(fitall)
separate <- with(H, ppm(P))
Fits <- hyperframe(Together=together, Separate=separate)
dr <- with(Fits, unlist(coef(Separate)) - unlist(coef(Together)))
dr
exp(dr)
@ 

One could also try deletion residuals, etc.

\subsection{Goodness-of-fit tests} 

\subsubsection{Quadrat count test}

The $\chi^2$ goodness-of-fit test based on quadrat counts is implemented
for objects of class \code{ppm} (in \code{quadrat.test.ppm})
and also for objects of class \code{mppm} (in \code{quadrat.test.mppm}).

This is a goodness-of-fit test for a fitted {\bf Poisson} point process
model only. The model could be uniform or non-uniform and the intensity
might depend on covariates. 

<<>>=
H <- hyperframe(X=waterstriders)

# Poisson with constant intensity for all patterns
fit1 <- mppm(X~1, H)
quadrat.test(fit1, nx=2)

# uniform Poisson with different intensity for each pattern
fit2 <- mppm(X ~ id, H)
quadrat.test(fit2, nx=2)
@ 

See the help for \code{quadrat.test.ppm} and \code{quadrat.test.mppm}
for further details.

\subsubsection{Kolmogorov-Smirnov test}

The Kolmogorov-Smirnov test of goodness-of-fit of a Poisson
point process model compares the observed and predicted
distributions of the values of a spatial covariate.

We want to test the null hypothesis $H_0$ that the observed point pattern
${\mathbf x}$ is a realisation from the Poisson process with intensity 
function $\lambda(u)$ (for locations $u$ in the window $W$).
Let $Z(u)$ be a given, real-valued covariate defined at each spatial location
$u$. Under $H_0$, the \emph{observed} values of $Z$ at the 
data points, $Z(x_i)$ for each $x_i \in {\mathbf x}$, are independent
random variables with common probability distribution function
\[
    F_0(z) = \frac{\int_W \lambda(u) \indicate{Z(u) \le z} \dee u}
                {\int_W \lambda(u)                       \dee u}.
\]
We can therefore apply the Kolmogorov-Smirnov test 
of goodness-of-fit. This compares the empirical cumulative distribution of
the observed values $Z(x_i)$ to the predicted c.d.f. $F_0$.

The test is implemented as \code{kstest.ppm}. The syntax is 

<<eval=FALSE>>=
kstest.mppm(model, covariate)
@ 

where \code{model} is a fitted model (of class \texttt{"mppm"})
and \code{covariate} is either 
\begin{itemize}
\item a \code{function(x,y)} making it possible to compute the value 
of the covariate at any location \code{(x,y)}
\item a pixel image containing the covariate values 
\item a list of functions, one for each row of the hyperframe of
original data
\item a list of pixel images, one for each row of the hyperframe of 
original data
\item a hyperframe with one column containing either functions or 
pixel images.
\end{itemize}

See Chapter 16 of \cite{baddrubaturn15} for further information.

\newpage

\addcontentsline{toc}{section}{Bibliography}

%\bibliography{%
%extra,%
%extra2,%
%biblio/badd,%
%biblio/bioscience,%
%biblio/censoring,%
%biblio/mcmc,%
%biblio/spatstat,%
%biblio/stat,%
%biblio/stochgeom%
%}

\begin{thebibliography}{1}

\bibitem{baddmollpake08}
A. Baddeley, J. M{\o}ller, and A.G. Pakes.
\newblock Properties of residuals for spatial point processes.
\newblock {\em Annals of the Institute of Statistical Mathematics},
  60:627--649, 2008.

\bibitem{TheBook}
A. Baddeley, E. Rubak, and R. Turner.
\newblock {\em Spatial Point Patterns: Methodology and Applications with R}.
\newblock Chapman \& Hall/CRC Press, 2015.

\bibitem{statpaper}
A. Baddeley, I. Sintorn, L. Bischof, R. Turner, and S. Heggarty.
\newblock Analysing designed experiments where the response is a spatial point
  pattern.
\newblock In preparation.

\bibitem{baddetal05}
A. Baddeley, R. Turner, J. M{\o}ller, and M. Hazelton.
\newblock Residual analysis for spatial point processes (with discussion).
\newblock {\em Journal of the Royal Statistical Society, series B},
  67(5):617--666, 2005.

\bibitem{chenetal08}
B.J. Chen, G.P. Leser, D. Jackson, and R.A. Lamb.
\newblock The influenza virus {M2} protein cytoplasmic tail interacts with the
  {M1} protein and influences virus assembly at the site of virus budding.
\newblock {\em Journal of Virology}, 82:10059--10070, 2008.

%#%^!ifdef RANDOMEFFECTS  
\bibitem{davigilt95}
M. Davidian and D.M. Giltinan.
\newblock {\em Nonlinear Mixed Effects Models for Repeated Measurement Data}.
\newblock Chapman and Hall, 1995.
%#%^!endif

\bibitem{digglangbene91}
P.J. Diggle, N. Lange, and F. M. Benes.
\newblock Analysis of variance for replicated spatial point patterns in
  clinical neuroanatomy.
\newblock {\em Journal of the {A}merican {S}tatistical {A}ssociation},
  86:618--625, 1991.

%#%^!ifdef RANDOMEFFECTS  
\bibitem{lairware82}
N.M. Laird and J.H. Ware.
\newblock Random-effects models for longitudinal data.
\newblock {\em Biometrics}, 38:963--974, 1982.
%#%^!endif

\bibitem{pent84}
A. Penttinen.
\newblock {\em Modelling Interaction in Spatial Point Patterns: Parameter
  Estimation by the Maximum Likelihood Method}.
\newblock Number 7 in {Jyv\"askyl\"a} Studies in Computer Science, Economics
  and Statistics. University of {Jyv\"askyl\"a}, 1984.

%#%^!ifdef RANDOMEFFECTS  
\bibitem{pinhbate00}
J.C. Pinheiro and D.M. Bates.
\newblock {\em Mixed-Effects Models in {S} and {S-PLUS}}.
\newblock Springer, 2000.
%#%^!endif

\end{thebibliography}

%\addcontentsline{toc}{section}{Index}

%\printindex

\end{document}