1

Good day,

I have been working through Baddeley et al. 2015 to fit a point process model to several point patterns using mppm {spatstat}. My point patterns are annual count data of large herbivores (i.e. point localities (x, y) of male/female animals * 5 years) in a protected area (owin). I have a number of spatial covariates e.g. distance to rivers (rivD) and vegetation productivity (NDVI).

Originally I fitted a model where herbivore response was a function of rivD + NDVI and allowed the coefficients to vary by sex (see mppm1 in reproducible example below). However, my annual point patterns are not independent between years in that there is a temporally increasing trend (i.e. there are exponentially more animals in year 1 compared to year 5). So I added year as a random effect, thinking that if I allowed the intercept to change per year I could account for this (see mppm2).

Now I'm wondering if this is the right way to go about it? If I was fitting a GAMM gamm {mgcv} I would add a temporal correlation structure e.g. correlation = corAR1(form=~year) but don't think this is possible in mppm (see mppm3)?

I would really appreciate any ideas on how to deal with this temporal correlation structure in a replicated point pattern with mppm {spatstat}.

Thank you very much

Sandra

# R version 3.3.1 (64-bit)
library(spatstat) # spatstat version 1.45-2.008

#### Simulate point patterns
# multitype Neyman-Scott process (each cluster is a multitype process)
nclust2 = function(x0, y0, radius, n, types=factor(c("male", "female"))) {
  X = runifdisc(n, radius, centre=c(x0, y0))
  M = sample(types, n, replace=TRUE)
  marks(X) = M
  return(X)
}

year1 = rNeymanScott(5,0.1,nclust2, radius=0.1, n=5)
# plot(year1)
#-------------------
year2 = rNeymanScott(10,0.1,nclust2, radius=0.1, n=5)
# plot(year2)
#-------------------
year2 = rNeymanScott(15,0.1,nclust2, radius=0.1, n=10)
# plot(year2)
#-------------------
year3 = rNeymanScott(20,0.1,nclust2, radius=0.1, n=10)
# plot(year3)
#-------------------
year4 = rNeymanScott(25,0.1,nclust2, radius=0.1, n=15)
# plot(year4)
#-------------------
year5 = rNeymanScott(30,0.1,nclust2, radius=0.1, n=15)
# plot(year5)

#### Simulate distance to rivers
line <- psp(runif(10), runif(10), runif(10), runif(10), window=owin())
# plot(line)
# plot(year1, add=TRUE)

#------------------------ UPDATE ------------------------#
#### Create hyperframe
#---> NDVI simulated with distmap to point patterns (not ideal but just to test)
hyp.years = hyperframe(year=factor(2010:2014),
                       ppp=list(year1,year2,year3,year4,year5),
                       NDVI=list(distmap(year5),distmap(year1),distmap(year2),distmap(year3),distmap(year4)),
                       rivD=distmap(line),
                       stringsAsFactors=TRUE)
hyp.years$numYear = with(hyp.years,as.numeric(year)-1)
hyp.years

#### Run mppm models
# mppm1 = mppm(ppp~(NDVI+rivD)/marks,data=hyp.years); summary(mppm1)
#..........................
# mppm2 = mppm(ppp~(NDVI+rivD)/marks,random = ~1|year,data=hyp.years); summary(mppm2)
#..........................
# correlation = corAR1(form=~year)
# mppm3 = mppm(ppp~(NDVI+rivD)/marks,correlation = corAR1(form=~year),use.gam = TRUE,data=hyp.years); summary(mppm3)

###---> Run mppm model with annual trend and random variation in growth
mppmCorr = mppm(ppp~(NDVI+rivD+numYear)/marks,random = ~1|year,data=hyp.years)
summary(mppm1)
Sands
  • 25
  • 5

1 Answers1

1

If there's a trend in population size over time, then it might make sense to include this trend in the systematic part of the model. I would suggest you add a new numeric variable NumYear to the data frame (eg giving the number of years since 2010). Then try adding simple trend terms such as +NumYear to the model formula (this would correspond to the exponential growth in population that you observed.) You can keep the 1|year random effect term which will then allow for random variation in population size around the long term growth trend.

There's no need to split the data patterns for each year into separate male and female patterns. The variable marks in the model formula can be used to specify any model that depends on sex.

I'm pretty sure that mppm with use.gam=TRUE does not recognise the argument correlation and this is probably just ignored. (It depends what happens inside gam).

Adrian Baddeley
  • 1,956
  • 1
  • 8
  • 7
  • Thank you Prof. Baddeley for your help. I have updated my script above as you suggested. I had split male and females because I could not get *marks* to work in the formula. I tried both *marks* and *'marks'* but I continue to get the following error: *invalid model formula in ExtractVars*. Any suggestions? – Sands Jun 27 '16 at 06:58
  • Thanks for spotting this. It must be a bug. Do you agree @adrian-baddeley? I have made a small correction in the `mppm` code which hopefully solves the problem with marks in the formula. You should use `marks` without any quotes. You can test this new version by installing a temporary branch of `spatstat` called "marked_mppm" (requires `devtools` package): `devtools::install_github("spatstat/spatstat", ref = "marked_mppm")` – Ege Rubak Jun 27 '16 at 07:47
  • Excellent! Thank you very much. I have updated the above script accordingly. – Sands Jun 27 '16 at 16:17
  • Dang. Yes, that's a bug. Ege's fix is correct. [The bug fix will be incorporated in the main branch of the development version of spatstat.] – Adrian Baddeley Jun 28 '16 at 04:15