[R-br] Modelo de regressão envolvendo rasters e vetores numéricos

[Jefferson Ferreira-Ferreira]

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 em 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*
>
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20170809/6b919f5b/attachment.html>