Pessoal, acho que não formulaei bem o problema. Por favor, considerem a mensagem abaixo como a correta. Agradeço se alguém puder dar alguma indicação. Abraços. library(raster) library(logistf) # Criando 2 rasters para reproduzir o que estou enfrentando r <- raster(nrow=10, ncol=10) # Variável-raster resposta (binária) com 9 bandas/layers s1 <- lapply(1:9, function(i) setValues(r, sample(0:1,ncell(r),replace = T))) s1 <- stack(s1) # um raster-variável preditora val <- sample(0:60,ncell(r),replace = T) s2 <- raster(nrow=10, ncol=10,vals=val) plot(s1) plot(s2) # Uma segunda variável preditora armazenada como vetor numérico: 9 valores assim como existem 9 bandas/layres em s1 (não é coincidência) exp_2 <- c(27.00,30.02,31.07,32.72,33.73,35.12,35.65,36.06,38.32) Agora, eu quero ajustar o modelo usando o logistf (eu tenho motivos estatísticos para isso, não reproduzido aqui) usando a função calc do pacote raster. É aí que o mistério mora. O raciocínio na formulação do problema é que cada célula/pixel seguiria a segunite fórmula glm: s1 ~ 27.00 + pixel correspondente em s2 + 27.00:pixel correspondente em s2 s1 ~ 30.02 + pixel correspondente em s2 + 30.02:pixel correspondente em s2 s1 ~ 31.07 + pixel correspondente em s2 + 31.07:pixel correspondente em s2 ... e assim por diante para todos os 9 valores de exp_2. Eu tentei algo assim: fun <- function (x) { logistf (x ~ exp_2 + s2 + exp_2:s2)$coeficients } coefs <- calc(s1, fun) Mas estava claro que não funcionaria. Provavelmente por causa das diferentes dimensões/número de rasterlayers neste conjunto de dados. Estou recebendo a mensagem de erro: Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : cannot use this function *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* Em 9 de agosto de 2017 16:07, Jefferson Ferreira-Ferreira <jecogeo@gmail.com
escreveu:
Pessoal,
Estou tentando ajustar a Estimativa de Máxima Máxima Penalizada de Firth para a regressão logística implementada no pacote logistf para um conjunto de rasters (imagens), mas estou com problemas para implementá-lo usando rasters e vetores numéricos. Para ilustrar o problema:
library(raster) library(logistf)
# Criando 2 rasters para reproduzir o que estou enfrentando r <- raster(nrow=10, ncol=10)
# Variável-raster resposta (binária) com 9 bandas/layers s1 <- lapply(1:9, function(i) setValues(r, sample(0:1,ncell(r),replace = T))) s1 <- stack(s1)
# um raster-variável preditora val <- sample(0:60,ncell(r),replace = T) s2 <- raster(nrow=10, ncol=10,vals=val)
plot(s1) plot(s2)
# Uma segunda variável preditora armazenada como vetor numérico: 9 valores assim como existem 9 bandas/layres em s1 (não é coincidência) exp_2 <- c(27.00,30.02,31.07,32.72,33.73,35.12,35.65,36.06,38.32)
Agora, eu quero ajustar o modelo usando o logistf (eu tenho motivos estatísticos para isso, não reproduzido aqui) usando a função calc do pacote raster. É aí que o mistério mora.
O raciocínio na formulação do problema é que cada célula/pixel seguiria a segunite fórmula glm:
s1/layer1 ~ 27.00 + pixel correspondente em s2 + 27.00:pixel correspondente em s2 s1/layer2 ~ 30.02 + pixel correspondente em s2 + 30.02:pixel correspondente em s2 s1/layer3 ~ 31.07 + pixel correspondente em s2 + 31.07:pixel correspondente em s2
... e assim por diante para todas as 9 bandas/layers do meu raster-variável resposta s1, que são pareadas com valores de exp_2.
Eu tentei algo assim:
fun <- function (x) { logistf (x ~ exp_2 + s2 + exp_2:s2)$coeficients } coefs <- calc(s1, fun)
Mas estava claro que não funcionaria. A parte complicada é infomar ao script que eu quero que cada valor de exp_2 seja usado com cada camada raster de s1 para este modelo.
Qualquer idéia será muito útil. Idéias?
Obrigado e abraços.
*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*