2

I'm running a spatial model using spatstat. Following model fitting, I would like to compare predictions at a few distances and say whether they are significantly different from each other. In the basis of the dataset, I have a bunch of points around the 0, 0 origin. Each point has x, y coordinates, distance from the origin, and cardinal angle (as north / east / west / south).

Following the setup from the answer to this question:

library(spatstat)
library(ggplot2)

# making up some data
set.seed(42)
X <- runifdisc(2000)
plot(X)
W <- Window(X)

rad <- as.im(function(x,y){sqrt(x^2+y^2)}, W)
ang <- as.im(atan2, W)
plot(solist(rad, ang), main = "")

# assigning directions
north <- ang < 45/180*pi & ang > -45/180*pi
east <- ang > 45/180*pi & ang < 135/180*pi
west <- ang < -45/180*pi & ang > -135/180*pi
south <- ang < -135/180*pi | ang > 135/180*pi
# create and run a model
lam <- 2000*exp(-2*rad - rad*north - 3*rad*west)
set.seed(42)
X2 <- rpoispp(lam)[W]

mod2 <- ppm(X2 ~ rad*west + rad*south +rad*east)
plot(predict(mod2))
plot(X2, add = TRUE, col = rgb(.9,.9,.9,.5))

Predict on a grid for north and west points. 

preds1 <- data.frame(Angle = "North", x = 0, y = seq(0.1, 1, 0.1))
preds2 <- data.frame(Angle = "West", y = 0, x = seq(-1, -0.1, 0.1)) 

preds <- rbind(preds1, preds2)
preds$Pred <- predict(mod2, locations = preds)
temp <- predict(mod2, locations = preds, interval = "confidence")
preds$Lower <- temp[1:nrow(preds)]
preds$Upper <- temp[(nrow(preds) + 1) : length(temp)]
preds$Distance <- ifelse(preds$x == 0, abs(preds$y), abs(preds$x))

ggplot(preds) +
   geom_ribbon(aes(x = Distance, ymin = Lower, ymax = Upper, group = Angle), alpha = 0.1) +
   geom_line(aes(x = Distance, y = Pred, colour = Angle), size = 0.75)

As the next step, I want to be able to say that north and west values are significantly different at distances x, y, z but not a, b, c... How do I get there?

user2602640
  • 640
  • 6
  • 21
  • Do I understand the question correctly: you want to say that they are different in the approx. range [0.12,1], but not in [0,0.12)? I.e., the goal is to determine where the confidence intervals do and don't overlap? – Julius Vainora Feb 04 '19 at 12:04
  • Yup, that looks right. – user2602640 Feb 04 '19 at 16:38
  • Looking at the defined distances (0.1, 0.2, ...), only 0.1 would be not significantly different. Is that all you want? Or do you want to *interpolate* the results as to actually get an interval as in my previous comment, [0,12]? – Julius Vainora Feb 04 '19 at 19:52
  • This example is a toy dataset - I'm looking for a general solution I can apply to my own data. I'd have to see how the interpolated results differ from the point-wise ones, can't quite picture it. In my mind, I'm looking at density, so point-wise comparisons make sense to me, but I'd be really curious to see both approaches. – user2602640 Feb 04 '19 at 20:07
  • Any luck? You got my hopes up :-) – user2602640 Feb 06 '19 at 10:35
  • The pointwise version is easy: for instance, we may first order the data with `preds <- preds[order(preds$Angle, preds$Distance), ]` and then use `preds$Distance[1:10][head(preds$Lower, 10) < tail(preds$Upper, 10)]`. I didn't write an answer as the interpolated version is going to be more involved.. The first step is to find where lower/upper limits intersect. The issue is that, I guess, North isn't always above West. So, that would require more care. In your example, we can use `rootSolve` to find multiple roots of a function f=f1-f2, where... – Julius Vainora Feb 06 '19 at 10:57
  • ...`f1 <- approxfun(preds$Distance[1:10], preds$Lower[1:10])` and `f2 <- approxfun(preds$Distance[11:20], preds$Upper[11:20])`. Then `multiroot(function(x) fun1(x) - fun2(x), start = 0.1)` finds 0.1909556 as the root and the answer would be [0.1909]. – Julius Vainora Feb 06 '19 at 10:58
  • I suggest to explain clearly in your answer what you want as to help people to avoid wasting time. And to think whether that's the right place to ask your question: is it methodological or about programming? – Julius Vainora Feb 09 '19 at 13:05
  • The second sentence in my question was "I would like to compare predictions at a few distances and say whether they are significantly different from each other". Which is what I want. I'm sorry it ended up being a waste of your time. I do think it's a programming question - I want a solution similar to `emmeans`, but for spatial models, at which point it becomes an application issue. – user2602640 Feb 09 '19 at 14:28

0 Answers0