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 R-sig-Geo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo