ppm                 package:spatstat                 R Documentation

_F_i_t _P_o_i_n_t _P_r_o_c_e_s_s _M_o_d_e_l _t_o _D_a_t_a

_D_e_s_c_r_i_p_t_i_o_n:

     Fits a point process model to an observed point pattern

_U_s_a_g_e:

      ppm(Q, trend=~1, interaction=NULL, covariates=NULL,
      correction="border", rbord=0, use.gam=FALSE, method="mpl",
      forcefit=FALSE, nsim=100, nrmh=1e5,
             start=list(n.start=X$n),
              control=list(nrep=nrmh),
              verb=TRUE)

_A_r_g_u_m_e_n_t_s:

       Q: A data point pattern (of class '"ppp"') to which the model
          will be fitted, or a quadrature scheme (of class '"quad"')
          containing this pattern. 

   trend: An R formula object specifying the spatial trend to be
          fitted.  The default formula, '~1', indicates the model is
          stationary and no trend is to be fitted.  

interaction: An object of class '"interact"' describing the point
          process interaction structure, or 'NULL' indicating that a
          Poisson process (stationary or nonstationary) should be
          fitted. 

covariates: The values of any spatial covariates (other than the
          Cartesian coordinates) required by the model. Either a data
          frame, or a list of images. See Details. 

correction: The name of the edge correction to be used. The default  is
          '"border"' indicating the border correction. Other
          possibilities may include '"Ripley"', '"isotropic"',
          '"translate"' and '"none"', depending on the  'interaction'. 

   rbord: If 'correction = "border"' this argument specifies the
          distance by which the window should be eroded for the border
          correction. 

 use.gam: Logical flag; if 'TRUE' then computations are performed using
          'gam' instead of 'glm'. 

  method: The method used to fit the model. Options are  '"mpl"' for
          the method of Maximum PseudoLikelihood, and '"ho"' for the
          Huang-Ogata approximate maximum likelihood method. 

forcefit: Logical flag for internal use. If 'forcefit=FALSE', some
          trivial models will be fitted by a shortcut. If
          'forcefit=TRUE', the generic fitting method will always be
          used.  

    nsim: Number of simulated realisations to generate (for
          'method="ho"') 

    nrmh: Number of Metropolis-Hastings iterations for each simulated
          realisation (for 'method="ho"') 

start,control: Arguments passed to 'rmh' controlling the behaviour of
          the Metropolis-Hastings algorithm (for 'method="ho"') 

    verb: Logical flag indicating whether to print progress reports
          (for 'method="ho"') 

_D_e_t_a_i_l_s:

     This function fits a point process model to an observed point
     pattern. The model may include spatial trend, interpoint
     interaction, and dependence on covariates.

     _b_a_s_i_c _u_s_e: In basic use, 'Q' is a point pattern dataset (an object
          of class '"ppp"') to which we wish to fit a model.

          The syntax of 'ppm()' is closely analogous to the R functions
          'glm' and 'gam'. The analogy is:

            *glm*      *ppm*
            'formula'  'trend'
            'family'   'interaction'

          The point process model to be fitted is specified by the 
          arguments 'trend' and 'interaction' which are respectively
          analogous to the 'formula' and 'family' arguments of glm(). 

          Systematic effects (spatial trend and/or dependence on 
          spatial covariates) are specified by the argument 'trend'.
          This is an R formula object, which may be expressed in terms
          of the Cartesian coordinates 'x', 'y', the marks 'marks', or
          the variables in 'covariates' (if supplied), or both. It
          specifies the *logarithm* of the first order potential of the
          process. The formula should not use any names beginning with
          '.mpl' as these are reserved for internal use. If 'trend' is
          absent or equal to the default, '~1', then the model to be
          fitted is stationary (or at least, its first order  potential
          is constant). 

          Stochastic interactions between random points of the point
          process are defined by the argument 'interaction'. This is an
          object of class '"interact"' which is initialised in a very
          similar way to the usage of family objects in 'glm' and
          'gam'. The families currently available are: 'Poisson',
          'Strauss', 'StraussHard', 'MultiStrauss', 'MultiStraussHard',
          'Softcore', 'DiggleGratton', 'Pairwise', 'PairPiece',
          'Geyer', 'LennardJones', 'Saturated', 'OrdThresh', and 'Ord'.
          See the examples below.

          If 'interaction' is missing or 'NULL', then the model to be
          fitted has no interpoint interactions, that is, it is a
          Poisson process (stationary or nonstationary according to
          'trend'). In this case the method of maximum pseudolikelihood
          coincides with maximum likelihood. 

          The fitted point process model returned by this function can
          be printed  (by the print method 'print.ppm') to inspect the
          fitted parameter values. If a nonparametric spatial trend was
          fitted, this can be extracted using the predict method
          'predict.ppm'.

     _M_o_d_e_l_s _w_i_t_h _c_o_v_a_r_i_a_t_e_s: To fit a model involving spatial
          covariates other than the Cartesian coordinates x and y, the
          values of the covariates should be supplied in the argument
          'covariates'.  Note that it is not sufficient to have
          observed the covariate only at the points of the data point
          pattern;  the covariate must also have been observed at other
           locations in the window.

          The argument 'covariates', if supplied, should be either a
          list of images, or a data frame.

          Typically 'covariates' is a list of images, each image giving
          the values of a spatial covariate at  a fine grid of
          locations.  The argument 'covariates' should then be a named
          list, with the names corresponding to the names of covariates
          in the model formula 'trend'. Each entry in the list must be
          an image object (of class '"im"', see 'im.object'). The
          software will look up the pixel values of each image at the
          required locations (quadrature points). In this case the
          argument 'Q' may be either a quadrature scheme or a point
          pattern.

          The variable names 'x', 'y' and 'marks' are reserved for the
          Cartesian  coordinates and the mark values, and these should
          not be used for variables in 'covariates'.

          If 'covariates' is a data frame, 'Q' must be a quadrature
          scheme (see under Quadrature Schemes below). Then
          'covariates' must have as many rows as there are points in
          'Q'. The ith row of 'covariates' should contain the values of
          spatial variables which have been observed at the ith point
          of 'Q'. 

     _Q_u_a_d_r_a_t_u_r_e _s_c_h_e_m_e_s: In advanced use, 'Q' may be a `quadrature
          scheme'. This was originally just a technicality but it has
          turned out to have practical uses, as we explain below.

          Quadrature schemes are required for our implementation of the
          method of maximum pseudolikelihood.  The definition of the
          pseudolikelihood involves an integral over the spatial window
          containing the data. In practice this integral must be
          approximated by a finite sum over a set of quadrature points.
          We use the technique of Baddeley and Turner (2000), a
          generalisation of the Berman-Turner (1992) device. In this
          technique the quadrature points for the numerical
          approximation include all the data points (points of the
          observed point pattern) as well as additional `dummy' points. 

          A quadrature scheme is an object of class '"quad"' (see
          'quad.object') which specifies both the data point pattern
          and the dummy points for the quadrature scheme, as well as
          the quadrature weights associated with these points. If 'Q'
          is simply a point pattern (of class '"ppp"', see
          'ppp.object') then it is interpreted as specifying the data
          points only; a set of dummy points specified by
          'default.dummy()' is added, and the default weighting rule is
          invoked to compute the quadrature weights.

          Finer quadrature schemes (i.e. those with more dummy points)
          generally yield a better approximation, at the expense of
          higher computational load.  Complete control over the
          quadrature scheme is possible. See 'quadscheme' for an
          overview.

          A practical advantage of quadrature schemes arises when we
          want to fit a model involving covariates (e.g. soil pH).
          Suppose we have only been able to observe the covariates at a
          small number of locations. Suppose 'cov.dat' is a data frame
          containing the values of the covariates at the data points
          (i.e. 'cov.dat[i,]' contains the observations for the 'i'th
          data point) and 'cov.dum' is another data frame (with the
          same columns as 'cov.dat') containing the covariate values at
          another set of points whose locations are given by the point
          pattern 'Y'. Then setting 'Q = quadscheme(X,Y)' combines the
          data points and dummy points into a quadrature scheme, and 
          'covariates = rbind(cov.dat, cov.dum)' combines the covariate
          data frames. We can then fit the model by calling 'ppm(Q,
          ..., covariates)'.

     _M_o_d_e_l-_f_i_t_t_i_n_g _t_e_c_h_n_i_q_u_e: The model may be fitted either by the
          method of maximum pseudolikelihood (Besag, 1975) or by the
          approximate maximum likelihood method of Huang and Ogata
          (1999). Maximum pseudolikelihood is much faster, but has
          poorer statistical properties. 

          In either case, the algorithm will begin by fitting the model
          by maximum pseudolikelihood. By default the algorithm returns
          the maximum pseudolikelihood fit.

          Maximum pseudolikelihood is equivalent to maximum likelihood
          for Poisson point processes. 

          Note that the method of maximum pseudolikelihood is believed
          to be inefficient and biased for point processes with strong
          interpoint interactions. In such cases, the Huang-Ogata
          approximate maximum likelihood method should be used,
          although maximum pseudolikelihood may also be used profitably
          for model selection in the initial phases of modelling.

     _H_u_a_n_g-_O_g_a_t_a _m_e_t_h_o_d: If 'method="ho"' then the model will be fitted
          using the Huang-Ogata (1999) approximate maximum likelihood
          method. First the model is fitted by maximum pseudolikelihood
          as described above, yielding an initial estimate of the
          parameter vector theta0. From this initial model, 'nsim'
          simulated realisations are generated. The score and Fisher
          information of the model at theta=theta0 are estimated from
          the simulated realisations. Then one step of the Fisher
          scoring algorithm is taken, yielding an updated estimate
          theta1. The corresponding model is returned.

          Simulated realisations are generated using 'rmh'. The
          iterative behaviour of the Metropolis-Hastings algorithm is
          controlled by the arguments 'start' and 'control' which are
          passed to 'rmh'. As a shortcut, the argument 'nrmh'
          determines the number of Metropolis-Hastings iterations run
          to produce one simulated realisation (if 'control' is
          absent).

     _E_d_g_e _c_o_r_r_e_c_t_i_o_n Edge correction should be applied to the
          sufficient statistics of the model, to reduce bias. The
          argument 'correction' is the name of an edge correction
          method. The default 'correction="border"' specifies the
          border correction, in which the quadrature window (the domain
          of integration of the  pseudolikelihood) is obtained by
          trimming off a margin of width 'rbord' from the observation
          window of the data pattern. Not all edge corrections are
          implemented (or implementable) for arbitrary windows. Other
          options depend on the argument 'interaction', but these
          generally include '"periodic"' (the periodic or toroidal edge
          correction in which opposite edges of a rectangular window
          are identified) and '"translate"' (the translation
          correction, see Baddeley 1998 and Baddeley and Turner 2000).
          For pairwise interaction there is also Ripley's isotropic
          correction, identified by '"isotropic"' or '"Ripley"'.

_V_a_l_u_e:

     An object of class '"ppm"' describing a fitted point process
     model.

     The fitted parameters can be obtained just by printing this
     object. Fitted spatial trends can be extracted using the 'predict'
     method for this object (see 'predict.ppm').

     See 'ppm.object' for details of the format of this object.

_W_a_r_n_i_n_g_s:

     The implementation of the Huang-Ogata method is experimental.

     See the comments above about the possible inefficiency and bias of
     the maximum pseudolikelihood estimator.

     The accuracy of the Berman-Turner approximation to the
     pseudolikelihood depends on the number of dummy points used in the
     quadrature scheme. The number of dummy points should  at least
     equal the number of data points.

     The parameter values of the fitted model do not necessarily
     determine a valid point process. Some of the point process models
     are only defined when the parameter values lie in a certain
     subset. For example the Strauss process only  exists when the
     interaction parameter gamma is less than or equal to 1,
     corresponding to a value of 'ppm()$theta[2]' less than or equal to
     '0'. The current version of 'ppm' maximises the pseudolikelihood
     without constraining the parameters, and does not apply any checks
     for sanity after fitting the model. 

     The 'trend' formula should not use any variable names beginning
     with the prefixes '.mpl' or 'Interaction' as these names are
     reserved for internal use. The data frame 'covariates' should have
     as many rows as there are points in 'Q'. It should not contain
     variables called 'x', 'y' or 'marks' as these names are reserved
     for the Cartesian coordinates and the marks.

     If the model formula involves one of the functions 'poly()',
     'bs()' or 'ns()' (e.g. applied to spatial coordinates 'x' and
     'y'), the fitted coefficients can be misleading. The resulting fit
     is not to the raw spatial variates ('x', 'x^2', 'x*y', etc.)  but
     to a transformation of these variates.  The transformation is
     implemented by 'poly()' in order to achieve better numerical
     stability. However the resulting coefficients are appropriate for
     use with the transformed variates, not with the raw variates.  
     This affects the interpretation of the constant term in the fitted
     model, 'logbeta'.  Conventionally, beta is the background
     intensity, i.e. the   value taken by the conditional intensity
     function when all predictors (including spatial or ``trend''
     predictors) are set equal to 0. However the coefficient actually
     produced is the value that the log conditional intensity takes
     when all the predictors,  including the _transformed_ spatial
     predictors, are set equal to '0', which is not the same thing.

     Worse still, the result of 'predict.ppm' can be completely wrong
     if the trend formula contains one of the functions 'poly()',
     'bs()' or 'ns()'. This is a weakness of the underlying function
     'predict.glm'. 

     If you wish to fit a polynomial trend,  we offer an alternative to
     'poly()', namely 'polynom()', which avoids the difficulty induced
     by transformations.  It is completely analogous to 'poly' except
     that it does not orthonormalise. The resulting coefficient
     estimates then have their natural interpretation and can be
     predicted correctly.  Numerical stability may be compromised.

     Values of the maximised pseudolikelihood are not comparable if
     they have been obtained with different values of 'rbord'.

_A_u_t_h_o_r(_s):

     Adrian Baddeley adrian@maths.uwa.edu.au <URL:
     http://www.maths.uwa.edu.au/~adrian/> and Rolf Turner
     rolf@math.unb.ca <URL: http://www.math.unb.ca/~rolf>

_R_e_f_e_r_e_n_c_e_s:

     Baddeley, A. and Turner, R. Practical maximum pseudolikelihood for
     spatial point patterns. _Australian and New Zealand Journal of
     Statistics_ *42* (2000) 283-322.

     Berman, M. and Turner, T.R.  Approximating point process
     likelihoods with GLIM. _Applied Statistics_ *41* (1992) 31-38.

     Besag, J. Statistical analysis of non-lattice data. _The
     Statistician_ *24* (1975) 179-195.

     Diggle, P.J., Fiksel, T., Grabarnik, P., Ogata, Y., Stoyan, D. and
     Tanemura, M. On parameter estimation for pairwise interaction
     processes. _International Statistical Review_ *62* (1994) 99-117.

     Huang, F. and Ogata, Y. Improvements of the maximum
     pseudo-likelihood estimators in various spatial statistical
     models. _Journal of Computational and Graphical Statistics_ *8*
     (1999) 510-530.

     Jensen, J.L. and Moeller, M. Pseudolikelihood for exponential
     family models of spatial point processes. _Annals of Applied
     Probability_ *1* (1991) 445-461.

     Jensen, J.L. and Kuensch, H.R.  On asymptotic normality of pseudo
     likelihood estimates for pairwise interaction processes, _Annals
     of the Institute of Statistical Mathematics_ *46* (1994) 475-486.

_S_e_e _A_l_s_o:

     'ppp',  'quadscheme',  'ppm.object', 'Poisson', 'Strauss',
     'StraussHard', 'MultiStrauss', 'MultiStraussHard', 'Softcore',
     'DiggleGratton', 'Pairwise', 'PairPiece', 'Geyer', 'LennardJones',
     'Saturated', 'OrdThresh', 'Ord'

_E_x_a_m_p_l_e_s:

      data(nztrees)
      ppm(nztrees)
      # fit the stationary Poisson process
      # to point pattern 'nztrees'

      Q <- quadscheme(nztrees) 
      ppm(Q) 
      # equivalent.

      ppm(nztrees, ~ x)
      # fit the nonstationary Poisson process 
      # with intensity function lambda(x,y) = exp(a + bx)
      # where x,y are the Cartesian coordinates
      # and a,b are parameters to be estimated

      ppm(nztrees, ~ polynom(x,2))
      # fit the nonstationary Poisson process 
      # with intensity function lambda(x,y) = exp(a + bx + cx^2)

      library(splines)
      ppm(nztrees, ~ bs(x,df=3))
      #       WARNING: do not use predict.ppm() on this result
      # Fits the nonstationary Poisson process 
      # with intensity function lambda(x,y) = exp(B(x))
      # where B is a B-spline with df = 3
      
      ppm(nztrees, ~1, Strauss(r=10), rbord=10)
      # Fit the stationary Strauss process with interaction range r=10
      # using the border method with margin rbord=10
       
      ppm(nztrees, ~ x, Strauss(13), correction="periodic")
      # Fit the nonstationary Strauss process with interaction range r=13
      # and exp(first order potential) =  activity = beta(x,y) = exp(a+bx)
      # using the periodic correction.

     # Huang-Ogata fit:
     ## Not run: ppm(nztrees, ~1, Strauss(r=10), rbord=10, method="mpl")

      # COVARIATES
      #
      X <- rpoispp(42)
      weirdfunction <- function(x,y){ 10 * x^2 + runif(length(x))}
      Zimage <- as.im(weirdfunction, unit.square())
      #
      # (a) covariate values in pixel image
      ppm(X, ~ y + Z, covariates=list(Z=Zimage))
      #
      # (b) covariate values in data frame
      Q <- quadscheme(X)
      xQ <- x.quad(Q)
      yQ <- y.quad(Q)
      Zvalues <- weirdfunction(xQ,yQ)
      ppm(Q, ~ y + Z, covariates=data.frame(Z=Zvalues))
      # Note Q not X

      ## MULTITYPE POINT PROCESSES ### 
      data(lansing)
      # Multitype point pattern --- trees marked by species


      # fit stationary marked Poisson process
      # with different intensity for each species
     ## Not run: ppm(lansing, ~ marks, Poisson())


      # fit nonstationary marked Poisson process
      # with different log-cubic trend for each species
     ## Not run: ppm(lansing, ~ marks * polynom(x,y,3), Poisson())

