[R-sig-Geo] Regression with rasters using calc::raster

[Jefferson Ferreira-Ferreira]

Dears,

I'm trying to fit the Firth's Penalized MLE GLM implemented in logistf
package to a set of rasters but I'm facing errors and problems I couldn't
realize until now.

# Lets generate 2 rasters reproducing what I'm facing
r <- raster(nrow=10, ncol=10)

# binary response raster-variable with 9 bands
s1 <- lapply(1:9, function(i) setValues(r, sample(0:1,ncell(r),replace =
T)))
s1 <- stack(s1)

# one explanatory raster-variable
val <- sample(0:60,ncell(r),replace = T)
s2 <- raster(nrow=10, ncol=10,vals=val)

plot(s1)
plot(s2)

# a second explanatory variable. Nine values
exp_2 <- c(27.00,30.02,31.07,32.72,33.73,35.12,35.65,36.06,38.32)

Now, I want to fit a model using Firth's Penalized MLE GLM implemented in
logistf (i have reasons for this not reproduced here) using calc from
raster package. That's where the mystery lives.

The rationale is each cell in:
s1/layer1 ~ 27.00 + corresponding cell in s2 + 27.00:corresponding cell in
s2
s1/layer2 ~ 30.02 + corresponding cell in s2 + 30.02:corresponding cell in
s2
s1/layer3 ~ 31.07 + corresponding cell in s2 + 31.07:corresponding cell in
s2
... and so on for all 9 bands of my response raster-variable, which are
paired with values from exp_2.

# I tried something like this:
fun <- function(x) { logistf(x ~ exp_2 + s2 + exp_2:s2)$coefficients }
coefs <- calc(s1,fun)

But it was clear it wouldn't work. The tricky part is to tell R I want each
value of exp_2 to be used with each rasterlayer of s1 for this model.

Any idea would be appreciated. Ideas?
Thanks in advance

*Jefferson Ferreira-Ferreira, **PhD (abd)*

*Geographer*



*Ecosystem Dynamics Observatory <http://tscanada.wix.com/ecodyn> -
EcoDyn/UNESP*
*Department of **Geography *
*Institute of Geosciences and Exact Sciences** (IGCE)*
*São Paulo State University (UNESP)*
*Rio Claro, SP - Brazil*

	[[alternative HTML version deleted]]