3

I am trying to do a regression analysis on a set of large rasters (254,004,000 cells each). Ultimately, I want to run something like the following (or a bit more complex, but let's start simple!):

model<-lm(dv ~ iv1+iv2+iv3... data=df,na.action=na.exclude)

where "dv" is the values from one raster and "iv1", "iv2" "iv3" ... are values from other rasters (up to 10 variables) with the same extent and resolution. It seems I should be able to do this out-of-memory using the Raster package, but I am confused how. Whether I create a brick, stack, or set of individual Raster objects, I cannot figure out how to send the variables to the lm function without using getValues and thus calling everything into memory (mine cannot even handle two variables).

A point in the right direction would be much appreciated!

Improve this question
5
  • I would take a moment and consider this in statistical terms. 1) you are effectively using the population and not a sample thus, negating the need for a regression. 2) using all the cells in a rasters is going to certainty add an unnecessary autocorrelation issue to a linear model. 3) In classical statistical terms, you will have a psuedoreplication (lack of independence) issue. 4) I highly doubt that you would meet iid assumptions. I would recommend taking a sample of the raster(s), use the sample data to build your regression model then estimate the model to your raster(s).
    – Jeffrey Evans
    Sep 26, 2013 at 23:26
  • Thanks, Jeffery, I appreciate the note. This would not be my final statistical product, but in conjunction with autocorrelation plots for each variable, I find it helps me with diagnostics. The answer below seems like it might be a fruitful path.
    – klwalker
    Sep 27, 2013 at 13:30
  • I beg to differ (slightly) with some of @Jeffrey Evans' points. First, regression for an entire population is meaningful: it describes relationships among variables. Second, autocorrelation is not necessarily a problem, but the advice to worry about it is excellent. There is a direct solution: tile your rasters. For each tile compute the mean, the count, and the [SSP matrix]. You can combine these statistics and proceed with the solution. There's no limit to the raster size this applies to. Another approach (using 2 rasters at a time) is given at stats.stackexchange.com/a/71257.
    – whuber
    Sep 27, 2013 at 15:56
  • I should be more specific. I do believe that regression approaches on rasters are useful in the context of "exploratory" analysis. One thought, have you considered an OLS rather than a straight linear model? The resulting residual error in OLS is a bit more robust to autocorrelation issues.
    – Jeffrey Evans
    Sep 27, 2013 at 19:22
  • Sorry, I should have been more specific as well. I do not necessarily want to pin myself to a linear model, I just thought that if I could get something to run with lm, I could carry it over to other similar packages (certainly not the most direct approach; a direct route to an OLS or other more robust method would be very welcome!) I have seen several examples of lm run in memory with the Raster package and get the impression that it can manage problems like these out of memory as well using a brick/stack object
    – klwalker
    Sep 27, 2013 at 20:18

2 Answers 2

Reset to default
2

The help for lm references biglm:

biglm in package biglm for an alternative way to fit linear models to large datasets.

The help pages for biglm indicate this package was developed for precisely such problems. The algorithm it references, AS274, is an updating procedure, allowing a solution based on a subset of the cases (cells) to be modified as additional cases are given.

Although this package appears to solve the problem, performing regression on such enormous datasets is (a) likely to be meaningless and (b) ignores opportunities to learn much more about the data. It is almost surely the case that posited relationships among the variables will change from one location to another. Why not capitalize on the size of the dataset, then, and conduct separate regressions within various windows or tiles of the data? For instance, you could tile your raster area into a 10 by 10 grid, reducing the size of each raster to less than three million cells, making in-memory calculations not only feasible but fast. If these regressions produce significantly different results you will have learned much (and will have avoided the error of combining them into one global regression); if they do not produce different results, you already have estimates of the global regression and you can justify regressing the entire dataset if you wish to do that.

It would likely be worthwhile to go even further and explore how the regressions vary with tile size. Much could be said about what tile sizes would be appropriate. I will limit my comments to just two simple ones. First, it would be a good idea to focus on sizes that are substantially larger than the longest range of spatial covariance of any of the variables. Second, the tiles need to be small enough to make computation practicable. There may be a wide range of choices between these two extremes.

Improve this answer
1

It looks like a Raster Stack object might be the most appropriate. I think this question might lead you in the right direction.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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