2

I am conducting a logistic regression analysis of spatial data for a project using STAN MCMC within R. The regression analysis itself is done, but now I need a way to apply the parameter estimates for different variables to raster data to generate prediction maps. Applying parameter estimates to continuous rasters in an ArcGIS Raster Calculator operation is straightforward (e.g., Intercept + 0.44*elevation - 0.23*dist_to_water...), but incorporating categorical rasters is more difficult. STAN requires categorical variables to be split up into a series of dummy variables, so my categorical rasters (e.g., native veg, surface geology, erosion class) need to be split up into a series of presence/absence (0/1) rasters for each value. For example, Native Vegetation has four categories, so I need four different 0/1 rasters that correspond to the distribution of each vegetation type. This is a similar problem to that discussed here.

However, I have six categorical rasters that each include a large number of values. As such, I would like to avoid having to manually create 0/1 rasters for each value using tools in ArcGIS like Extract by Attributes, Reclassify, or Conditional statements (e.g., Con(native_veg, 1, 0, "VALUE = 1")).

Is there a function in R that can iterate through unique raster values and create new 0/1 rasters for each?

1

The raster package has a function to do this for you in one line. Use layerize():

library(raster)
# make an example raster
r <- raster(nrow=100, ncol=100)
r[] <- round(runif(ncell(r),1,4),0)

# create presence/absence rasters, stored in a RasterBrick. 
r_dummy <- layerize(r)
plot(r_dummy[[1]])

Raster for "1" values

  • Very cool. This is a much simpler solution, so I accepted your answer instead of the original. Thanks! – lambertj Sep 17 at 14:18
3

Are you sure that the dummy coding in STAN requires separate binomial parameters for the estimate? For efficiency sake, the dummy coding for say, Native Vegetation should be a factor of [1,2,3,4] and not four separate parameters. If this is, in fact, the case it is trivial to do this in R.

library(raster)
r <- raster(nrow=100, ncol=100)
  r[] <- round(runif(ncell(r),1,4),0)

r1=r; r2=r; r3=r; r4=r # make copies of original to modify
  r1[] <- ifelse(r[] == 1, 1, 0)
  r2[] <- ifelse(r[] == 2, 1, 0)
  r3[] <- ifelse(r[] == 3, 1, 0)
  r4[] <- ifelse(r[] == 4, 1, 0)

plot(stack(r1,r2,r3,r4))

You can extend this into processing all of your rasters in a for loop.

First, as an example, create raster stack with 5 layers and different unique values ranging 1 - 5. This should emulate you problem since each raster has an different number of levels. This simulated data is, of course, representing a raster stack or brick object of you nominal covariates that need to be converted to binomial. Your raster data can be read into this object class using the raster::stack or raster::brick functions.

r <- stack(brick(array(runif(100 * 100 * 5), dim=c(100, 100, 5))))
  for(i in 1:nlayers(r)) { r[[i]][] <- round(runif(ncell(r), 
                           sample(1:2,1),sample(3:5,1) ),0) }
names(r) <- paste0("parameter",1:5)

Now we can define a double loop that looks at the unique values for each raster and then loops through them to create binary rasters for each rasters unique levels. The resulting object (given the below code) will be binary.rasters.

binary.rasters <- stack()
rnames <- vector() 
  for(i in 1:nlayers(r)) {
    for(j in 1:length(unique(r[[i]][]))) {
      u <- unique(r[[i]][])
      b <- r[[i]] 
      b[] <- ifelse(b[] == u[j], 1, 0)
      binary.rasters <- addLayer(binary.rasters, b)
      rnames <- append(rnames, paste(names(r)[i], paste("level",u[j],sep="-"), sep="_") )  
    }
  }
( names(binary.rasters) <- rnames )

I added a vector rnames that tracks the raster name and appends it with the level in each sub loop. This can then be added as the names of the resulting raster stack. This will allow one to know what raster and level a given binary raster is representing however, does not necessary correspond with the parameter names used in the original model. Although, you can easily rename the elements in this vector.

I have seen this type of dummy coding in occupancy modeling (eg., software Presence) and it sure eats up parameter space. One has to wonder if it would not be prudent to form explicit hypothesis around a specific level in a given categorical variable and not think of it as a single independent variable, because that is sure not how it is behaving in the model. In this case each level in the variable is truly a separate parameter. If there are only a level or two in a given independent nominal variable that you hypothesize is going to effect your process, why include all of the levels? By reducing the parameter space you are going to get much more relevant AIC values when evaluating competing models.

  • Unfortunately, STAN can only accept real or integer variables. All my categorical variables are factors in the original data, but I broke them into separate dummy variables using: library(cobalt) dcc.s.dummy <- splitfactor(dcc.s) If you try to feed unsplit factors into an rSTAN model, it returns the error: In FUN(X[[i]], ...) : data with name VARIABLE is not numeric and not used. – lambertj Jun 20 '18 at 20:36
  • I completely see your point about only using specific levels of categorical variables that test an explicit hypothesis. I've already done that for things like Native Vegetation and Surface Geology. Unfortunately, the structure of the data I have available contains a number of other variables that probably could (and should) be continuous, but are only collected as categorical (e.g., erosion class instead of proportion of topsoil eroded). There isn't much I can do about that-- that's just the data I have to work with. – lambertj Jun 20 '18 at 20:36
  • I modified your code above to test on one of my categorical raster layers (glacial drift thickness): drift_binary <- stack() for(j in 1:length(unique(drift_thick))) { u <- unique(drift_thick[]) b <- drift_thick b[] <- ifelse(b[] == u[j], 1, 0) drift_binary <- addLayer(drift_binary, b) } writeRaster(drift_binary, 'drift.tif', bylayer=TRUE) – lambertj Jun 20 '18 at 22:18
  • This does work and produces a series of binary 0/1 rasters. However, the names of the layers in the raster stack (drift_thick.1.1, drift_thick.2.1, etc.) and output rasters (drift_1, drift_2, etc.) don't match the actual parent variable levels. Is there a way to force the layers in the raster stack to take on the names of the variable levels? – lambertj Jun 20 '18 at 22:18
  • In the addLayer function you could try assign to create a different variable name or could create an empty vector and add to it using paste to create a vector of names that is then assigned back to the final object. I will append my answer when i am back at a computer. – Jeffrey Evans Jun 20 '18 at 22:49

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.