rmhmodel              package:spatstat              R Documentation

_D_e_f_i_n_e _P_o_i_n_t _P_r_o_c_e_s_s _M_o_d_e_l _f_o_r _M_e_t_r_o_p_o_l_i_s-_H_a_s_t_i_n_g_s _S_i_m_u_l_a_t_i_o_n.

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

     Builds a description of a point process model for use in
     simulating the model by the Metropolis-Hastings algorithm.

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

        rmhmodel(model)
        rmhmodel(..., cif=NULL, par=NULL, w=NULL, trend=NULL, types=NULL)

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

   model: An existing description of the model in some other format.
          Incompatible with the arguments listed below. 

     ...: There should be no other arguments.

     cif: Character string specifying the choice of model

     par: Parameters of the model

       w: Spatial window in which to simulate

   trend: Specification of the trend in the model

   types: A vector of factor levels defining the possible marks, for a
          multitype process. 

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

     Simulated realisations of many point process models can be
     generated using the Metropolis-Hastings algorithm 'rmh'. The
     algorithm requires the model to be specified in a particular
     format.

     This function 'rmhmodel' creates a description of the point
     process model in the form required by 'rmh', and checks that the
     parameters of the model are valid.

     The point process model should be specified *either* by the first
     argument 'model' *or* by the other arguments 'cif', 'par' etc.

     If 'model' is a fitted point process model (object of class
     '"ppm"' obtained by a call to the model-fitting function 'ppm')
     then a description of the point process model will be extracted
     from this object.

     If 'model' is a list, then it should have components named  'cif',
     'par' and optionally 'w', 'trend', 'types' with the same
     interpretation as described below.

     The argument 'cif' is a character string specifying the choice of
     interpoint interaction for the point process. The current options
     are

     ''_s_t_r_a_u_s_s'' The Strauss process

     ''_s_t_r_a_u_s_h'' The Strauss process with hard core

     ''_s_f_t_c_r'' The Softcore process

     ''_s_t_r_a_u_s_s_m'' The multitype Strauss process

     ''_s_t_r_a_u_s_h_m'' Multitype Strauss process with hard core

     ''_d_g_s'' Diggle, Gates and Stibbard (1987) process

     ''_d_i_g_g_r_a'' Diggle and Gratton (1984) process

     ''_g_e_y_e_r'' Saturation process (Geyer, 1999).

     ''_l_o_o_k_u_p'' General isotropic pairwise interaction process, with
          the interaction function specified via a ``lookup table''.

     The argument 'par' supplies parameter values appropriate to the
     conditional intensity function being invoked. These are:

     _s_t_r_a_u_s_s: (Strauss process.) A *named* vector with components
          'beta,gamma,r' which are respectively the ``base'' intensity,
          the pair-wise interaction parameter and the interaction
          radius.  Note that 'gamma' must be less than or equal to 1.

     _s_t_r_a_u_s_h: (Strauss process with hardcore.) A *named* vector with
          entries 'beta,gamma,r,hc' where 'beta', 'gamma', and 'r' are
          as for the Strauss process, and 'hc' is the hardcore radius. 
          Of course 'hc' must be less than 'r'.

     _s_f_t_c_r: (Softcore process.) A *named* vector with components
          'beta,sigma,kappa'.  Again 'beta' is a ``base'' intensity.
          The pairwise interaction between two points u != v is

                      -(sigma/||u-v||)^(2/kappa)

          Note that it is necessary that 0 < kappa <1.

     _s_t_r_a_u_s_s_m: (Multitype Strauss process.) A *named* list with
          components

             *  'beta':  A vector of ``base'' intensities, one for each
                possible type.

             *  'gamma': A *symmetric* matrix of interaction
                parameters, with gamma_ij pertaining to the interaction
                between type i and type j.

             *  'radii': A *symmetric* matrix of interaction radii,
                with entries r_ij pertaining to the interaction between
                type i and type j.


     _s_t_r_a_u_s_h_m: (Multitype Strauss process with hardcore.) A *named*
          list with components 'beta' and 'gamma' as for 'straussm' and
          *two* ``radii'' components:

             *  'iradii': the interaction radii

             *  'hradii': the hardcore radii

          which are both symmetric matrices of nonnegative numbers. The
          entries of 'hradii' must be less than the corresponding
          entries of 'iradii'.

     _d_g_s: (Diggle, Gates, and Stibbard process. See Diggle, Gates, and
          Stibbard (1987)) A *named* vector with components 'beta' and
          'rho'.  This process has pairwise interaction function equal
          to

                   e(t) = sin^2((pi * t)/(2 * rho))

          for t < rho, and equal to 1 for t >= rho.

     _d_i_g_g_r_a: (Diggle-Gratton process. See Diggle and Gratton (1984) and
          Diggle, Gates and Stibbard (1987).) A *named* vector with
          components 'beta', 'kappa', 'delta' and 'rho'.  This process
          has pairwise interaction function e(t) equal to 0 for t <
          delta, equal to

                    ((t-delta)/(rho-delta))^kappa

          for delta <= t < rho, and equal to 1 for delta >= rho. Note
          that here we use the symbol kappa where Diggle, Gates, and
          Stibbard use beta since we reserve the symbol beta for an
          intensity parameter.

     _g_e_y_e_r: (Geyer's saturation process. See Geyer (1999).) A *named*
          vector with components 'beta', 'gamma', 'r', and 'sat'. The
          components 'beta', 'gamma', 'r' are as for the Strauss model,
          and 'sat' is the ``saturation'' parameter.  The model is
          Geyer's ``saturation'' point process model, a modification of
          the Strauss process in which we effectively impose an upper
          limit ('sat') on the number of neighbours which will be
          counted as close to a given point.

          Explicitly, a saturation point process with interaction
          radius r, saturation threshold s, and parameters beta and
          gamma, is the point process in which each point x[i] in the
          pattern X contributes a factor

                     beta gamma^min(s,t(x[i],X))

          to the probability density of the point pattern, where
          t(x[i],X) denotes the number of ``r-close neighbours'' of
          x[i] in the pattern X.


     _l_o_o_k_u_p: (Arbitrary pairwise interaction process with isotropic
          interaction.) A *named* list with components 'beta', 'r', and
          'h', or just with components 'beta' and 'h'.

          This model is the pairwise interaction process with an
          isotropic interaction given by any chosen function H. Each
          pair of points x[i], x[j] in the point pattern contributes a
          factor H(d(x[i],x[j])) to the probability density, where d
          denotes distance and H is the pair interaction function.

          The component 'beta' is a (positive) scalar which determines
          the ``base'' intensity of the process.

          In this implementation, H must be a step function. It is
          specified by the user in one of two ways.

             *  *as a vector of values:* If 'r' is present, then 'r' is
                assumed to  give the locations of jumps in the function
                H, while the vector 'h' gives the corresponding values
                of the function.

                Specifically, the interaction function H(t) takes the
                value 'h[1]' for distances t in the interval  '[0,
                r[1])'; takes the value 'h[i]' for distances t in the
                interval  '[r[i-1], r[i])' where i = 2, ..., n; and
                takes the value 1 for t >= r[n]. Here n denotes the
                length of 'r'.

                The components 'r' and 'h' must be numeric vectors of
                equal length. The 'r' values must be strictly positive,
                and  sorted in increasing order.

                The entries of 'h' must be non-negative.  If any entry
                of 'h' is greater than 1, then the entry 'h[1]' must be
                0 (otherwise the specified process is non-existent).

                Greatest efficiency is achieved if the values of 'r'
                are equally spaced.

                [*Note:* The usage of 'r' and 'h' has _changed_ from
                the previous usage in 'spatstat' versions 1.4-7 to
                1.5-1, in which ascending order was not required, and
                in which the first entry of 'r' had to be 0.]

             *  *as a stepfun object:* If 'r' is absent, then 'h' must
                be an object of class '"stepfun"' specifying a step
                function. Such objects are created by 'stepfun'. 

                The stepfun object 'h' must be right-continuous (which
                is the default using 'stepfun'.)

                The values of the step function must all be
                nonnegative. The values must all be less than 1 unless
                the function is identically zero on some initial
                interval [0,r). The rightmost value (the value of
                'h(t)' for large 't') must be equal to 1.

                Greatest efficiency is achieved if the jumps (the
                ``knots'' of the step function) are equally spaced.


     The optional argument 'trend' determines the spatial trend in the
     model, if it has one. It should be a function or image (or a list
     of such, if the model is multitype) to provide the value of the
     trend at an arbitrary point.

     _t_r_e_n_d _g_i_v_e_n _a_s _a _f_u_n_c_t_i_o_n: A trend function may be a function of
          any number of arguments, but the first two must be the
          eqn{x,y} coordinates of a point.  Auxiliary arguments may be
          passed to the 'trend' function at the time of simulation, via
          the '...' argument to 'rmh'.

          The function *must* be *vectorized*. That is, it must be
          capable of accepting vector valued 'x' and 'y' arguments. 
          Put another way, it must be capable of calculating the trend
          value at a number of points, simultaneously, and should
          return the *vector* of corresponding trend values.

     _t_r_e_n_d _g_i_v_e_n _a_s _a_n _i_m_a_g_e: An image (see 'im.object') provides the
          trend values at a grid of points in the observation window
          and determines the trend value at other points as the value
          at the nearest grid point.

     Note that the trend or trends must be *non-negative*; no checking
     is done for this.

     The optional argument 'w' specifies the window in which the
     pattern is to be generated.  If specified, it must be in a form
     which can be coerced to an object of class 'owin' by 'as.owin'.

     The optional argument 'types' specifies the possible types in a
     multitype point process. If the model being simulated is
     multitype, and 'types' is not specified, then this vector defaults
     to '1:ntypes' where 'ntypes' is the number of types.

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

     An object of class '"rmhmodel"', which is essentially a list of
     parameter values for the model.

     There is a 'print' method for this class, which prints a sensible
     description of the model chosen.

_W_a_r_n_i_n_g_s _i_n _R_e_s_p_e_c_t _o_f "_l_o_o_k_u_p":

     The syntax of 'rmh.default' in respect of the 'lookup' cif has
     _changed_ from the previous release of 'spatstat' (versions 1.4-7
     to 1.5-1). Read the _Details_ carefully. In particular it is now
     required that the first entry of the 'r' component of 'par' be
     _strictly positive_.  (This is the opposite of what was required
     in the previous release, which was that this first entry had to be
     0.)

     It is also now required that the entries of 'r' be sorted into
     ascending order.  (In the previous release it was assumed that the
     entries of 'r' and 'h' were in corresponding order and the two
     vectors were sorted commensurately.  It was decided that this is
     dangerous sand unnecessary.)

     Note that if you specify the 'lookup' pairwise interaction
     function via 'stepfun()' the arguments 'x' and 'y' which are
     passed to 'stepfun()' are slightly different from 'r' and 'h': 
     'length(y)' is equal to '1+length(x)'; the final entry of 'y' must
     be equal to 1 - i.e. this value is explicitly supplied by the user
     rather than getting tacked on internally.

     The step function returned by 'stepfun()' must be right continuous
     (this is the default behaviour of 'stepfun()') otherwise an error
     is given.

_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:

     Diggle, P. J. (2003) _Statistical Analysis of Spatial Point
     Patterns_ (2nd ed.) Arnold, London.

     Diggle, P.J. and Gratton, R.J. (1984) Monte Carlo methods of
     inference for implicit statistical models. _Journal of the Royal
     Statistical Society, series B_ *46*, 193 - 212.

     Diggle, P.J., Gates, D.J., and Stibbard, A. (1987) A nonparametric
     estimator for pairwise-interaction point processes. Biometrika
     *74*, 763 - 770. _Scandinavian Journal of Statistics_ *21*,
     359-373.

     Geyer, C.J. (1999) Likelihood Inference for Spatial Point
     Processes. Chapter 3 in  O.E. Barndorff-Nielsen, W.S. Kendall and
     M.N.M. Van Lieshout (eds) _Stochastic Geometry: Likelihood and
     Computation_, Chapman and Hall / CRC,  Monographs on Statistics
     and Applied Probability, number 80. Pages 79-140.

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

     'rmh', 'rmhcontrol', 'rmhstart', 'ppm', 'Strauss', 'Softcore',
     'StraussHard', 'MultiStrauss', 'MultiStraussHard',
     'DiggleGratton', 'PairPiece'

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

        # Strauss process:
        mod01 <- rmhmodel(cif="strauss",par=c(beta=2,gamma=0.2,r=0.7),
                      w=c(0,10,0,10))
        # Equivalent to:
        a <- list(cif="strauss",par=c(beta=2,gamma=0.2,r=0.7),
                      w=c(0,10,0,10))
        mod01 <- rmhmodel(a)

        # Strauss with hardcore:
        mod04 <- list(cif="straush",par=c(beta=2,gamma=0.2,r=0.7,hc=0.3),
                     w=owin(c(0,10),c(0,5)))
        mod04 <- rmhmodel(mod04)

        # Soft core:
        par3 <- c(0.8,0.1,0.5)
        w    <- square(10)
        mod07 <- rmhmodel(cif="sftcr",
                          par=c(beta=0.8,sigma=0.1,kappa=0.5),
                          w=w)
        
        # Multitype Strauss:
        beta <- c(0.027,0.008)
        gmma <- matrix(c(0.43,0.98,0.98,0.36),2,2)
        r    <- matrix(c(45,45,45,45),2,2)
        mod08 <- rmhmodel(cif="straussm",
                          par=list(beta=beta,gamma=gmma,radii=r),
                          w=square(250))
        # specify types
        mod09 <- rmhmodel(cif="straussm",
                          par=list(beta=beta,gamma=gmma,radii=r),
                          w=square(250),
                          types=c("A", "B"))

        
        # Multitype Strauss hardcore with trends for each type:
        beta  <- c(0.27,0.08)
        ri    <- matrix(c(45,45,45,45),2,2)
        rhc  <- matrix(c(9.1,5.0,5.0,2.5),2,2)
        tr3   <- function(x,y){x <- x/250; y <- y/250;
                                exp((6*x + 5*y - 18*x^2 + 12*x*y - 9*y^2)/6)
                              }
                              # log quadratic trend
        tr4   <- function(x,y){x <- x/250; y <- y/250;
                              exp(-0.6*x+0.5*y)}
                             # log linear trend
        mod10 <- rmhmodel(cif="straushm",par=list(beta=beta,gamma=gmma,
                      iradii=ri,hradii=rhc),w=c(0,250,0,250),
                      trend=list(tr3,tr4))

        # Lookup (interaction function h_2 from page 76, Diggle (2003)):
           r <- seq(from=0,to=0.2,length=101)[-1] # Drop 0.
           h <- 20*(r-0.05)
           h[r<0.05] <- 0
           h[r>0.10] <- 1
           mod17 <- rmhmodel(cif="lookup",par=list(beta=4000,h=h,r=r),w=c(0,1,0,1))

