I am wanting to estimate the relationship between two rasters in R, but am having some trouble thinking through how I should do it. I am providing a toy example.

For example, suppose x (e.g. NDVI) is giving rise to y, which is a variable that sums up multiple species responses to x.

x = runif(100)
y = .25 + (2.5*x) + (-2.5 * I(x^2))
eps = rnorm(length(x),mean = 0, sd = .1)
y = as.vector(y + eps)

Here are the rasters from the data

y_rast = raster(nrow=10,ncol=10)
values(y_rast) = y

x_rast = y_rast
values(x_rast) = x

What is the best way to express the relationship between y and x? I was thinking of sampling, but would like to include some kind of estimation of significance (used loosely here).


This is basically the whole of statistics.

Given you simulated the data from a quadratic with a specified random noise, ideally you'd get that back as your "expression of the relationship between y and x", and your significance would be in the uncertainty in the parameters of the quadratic and the residual noise estimate, using standard likelihood-based techniques.

From real data you don't know what generated it, so you make a guess based on exploratory analysis or first principles given the science or a hypothesis you want to test. Looking at your x-y plot, I'd go "hmm looks quadratic" and fit a quadratic. My first assumption would also be uncorrelated independent errors of equal magnitude.

Once fitted there's "goodness of fit" methods to tell you where your model assumptions go wrong. You might discover your errors are spatially correlated, or are much larger in one area than another, in which case you have to reformulate your model with these discoveries and fit again.

So the short answer is "statistics". Propose a model for the mean relationship, (ooh, quadratic...), model for the error (ooh, uncorrelated independent errors), then fit, refine, until you can't disprove yourself any more. Goes for spatial data, longitudinal data, experimental data - anything. Its statistics.

Is that too broad? I don't want to say "fit a polynomial model using glm(y~Poly(x, 2)) if your real data doesn't look like this. No, you need to go through the statistical process as outlined above.

  • Thanks. That is essentially what I was looking for. I am still trying to think of the best way to bootstrap a p-value, but maybe I just compare the model fits.
    – user44796
    Jun 14 at 15:45
  • You can't have a p-value without a statistical hypothesis, and you can't have that without a probabilistic model, which is what model fits are.
    – Spacedman
    Jun 14 at 16:51
  • That's true, my issue is that, as sample size increases the p-value decreases. So you can "cheat" and just sample more points. What I am thinking is that I am going to randomly sample points X times, and captures the parameter estimates. If the estimate doesn't cross 0 at 95% CI, should be fine.
    – user44796
    Jun 14 at 20:26
  • I don't see what you original question has got to do with sampling. Why sample? Why not use all your data?
    – Spacedman
    Jun 14 at 20:28
  • I wouldn't want to use all the points because of the relationship between p-values and sample size, and the spatial autocorrelation that comes with using all the points. I wouldn't necessarily call out someone that does it, I just don't feel comfortable doing it.
    – user44796
    Jun 15 at 12:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.