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]]

_______________________________________________

R-sig-Geo mailing list

[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo