[R-br] erro na validação cruzada com pacote geoR

Sábado Outubro 26 10:55:55 BRST 2013

Senhores, bom dia!

Agradeço novamente ao Prof. Paulo pela atenção dispensada.

Hélio, quanto ao script que eu tinha ficado de olhar, vou colocar minhas
considerações abaixo.

Procurei não mudar muito, mas acabei alterando um pouco. Os comentários
estão no corpo do script.

Lembrando sempre que também estou aprendendo sobre o assunto, então pode
ter muita coisa errada ou mal interpretada aí. Mas acho que pode te ajudar.

Outra coisa é que não tenho muito ideia do que os dados são, como foram
obtidos e as condições da área de onde vieram, o que pra mim é fundamental
pra análise. Pelos indícios nos dados (localização, estrutura), vou
arriscar que são dados da produção de café ou citros e que o terreno está
numa pendente da direita pra esquerda, com maior produção na parte baixa. O
que chamou a atenção é o formato da área que é bem particular e não parece
muito 'natural'.

Mas como o foco principal é o uso do R/geoR, vou tentar me deter nisso.

Use por sua própria conta e risco! :D

### <BEGIN>

require(geoR) #;require(tcltk2)
setwd("C:/LAB/RGIS/trend"); dir()
par.ori <- par(no.readonly = TRUE)

url1 <- 'https://dl.dropboxusercontent.com/u/117618178/helio/dados.txt'
url2 <- 'https://dl.dropboxusercontent.com/u/117618178/helio/recorte.txt'

### Como criar uma pasta pública no Dropbox...
browseURL('https://www.dropbox.com/help/16/pt_BR')

### 1. Carregando dados
------------------------------------------------------
dados <- read.geodata(url1, coords.col=1:2, data.col=3, sep="", header=T)
dados$borders <- read.table(url2, header=T)
summary(dados); length(dados$data)
ls(dados);ls.str(dados)

### 2. Visualizando dados
----------------------------------------------------

### Visão geral
par(mfrow=c(1,1)); par()$mfrow
points(dados)
points(dados, pt.div="quart")
points(dados, pt.div="quint")
### tem pelos três sub-áreas bem marcadas...

### Verificando tendências com o plot.geodata
--------------------------------
### Antes de partir pros variogramas, use o plot.geodata para pré-avaliar
tendências
### Bom para esclarecer opções de trend, pois notei certa confusão
(cte,1st,2nd)
### Recomendo ler o item 'trend' na ajuda da função plot.geodata
?plot.geodata

plot(dados, lowess=T)
plot(dados, lowess=T, trend='cte') ### mesmo que anterior
### linhas de comando anteriores são equivalentes: trend='cte' é 'default'
da geoR!
### nos gráf. 2 e 3, plota-se as coord vs. dados. Observar as distorções na
linha 'lowess'

plot(dados, lowess=T, trend=~coords)
plot(dados, lowess=T, trend='1st')
plot(dados, lowess=T, trend=~dados$coords[,1]+dados$coords[,2])
### as três linhas anteriores são equivalentes
### o uso da polinomial de primeira ordem melhora o comportamento do
histograma
### mas a linha de lowess continua com comportamento anômalo

plot(dados, lowess=T, trend='2nd')
plot(dados, lowess=T, trend=~I(coords)+I(coords^2)+I(coords[,1]*coords[,2]))
### as duas linhas anteriores são equivalentes (polinomal de segunda ordem)!
### não aparenta agregar nenhuma melhoria

### Mais sobre modelos de tendência com trend.spatial
------------------------
### Pra melhorar o entendimento dos modelos, dá pra dar uma olhada nas
matrizes
?trend.spatial
head(unclass(trend.spatial('cte', dados)))   ### constante
head(unclass(trend.spatial(~coords, dados))) ### polinomal 1a ordem, usando
coords
head(unclass(trend.spatial('1st', dados)))   ### mesmo que anterior
head(unclass(trend.spatial('2nd', dados)))   ### polinomal 2a ordem, usando
coords
head(unclass(trend.spatial(~coords[,1], dados))) ### em função de X
head(unclass(trend.spatial(~coords[,2], dados))) ### em função de Y
head(unclass(trend.spatial(~coords+I(coords[,1]^2), dados))) ### coords +
X**2
head(unclass(trend.spatial(~I(coords)+I(coords^2), dados)))
head(unclass(trend.spatial(~I(coords)+I(coords^2)+I(coords[,1]*coords[,2]),
dados))) ### mesmo que 2nd
head(unclass(trend.spatial(~I(coords)+I(coords[,1]^2), dados)))
head(unclass(trend.spatial(~I(coords^2), dados)))
head(unclass(trend.spatial(~coords[,1]*coords[,2], dados)))
head(unclass(trend.spatial(~I(coords[,1]*coords[,2]), dados)))
head(unclass(trend.spatial(~coords[,1]+I(coords[,1]^2)+coords[,2], dados)))

### Separando o drift com base na matriz trend.spatial
-----------------------
trend.matrix <- unclass(trend.spatial(trend='1st', geodata=dados));
nrow(trend.matrix)
trend.lm     <- lm(dados$data ~ trend.matrix)
trend.res    <- lm(dados$data ~ trend.matrix)$residuals
trend.fit    <- lm(dados$data ~ trend.matrix)$fitted.values

### 'Empacotando' os resíduos em um
geodata-----------------------------------
dados.res <- dados
names(trend.res) <- NULL; trend.res
dados.res$data <- trend.res; dados.res$data
str(dados.res)
### Pode usar dados.res pra fazer a análise variográfica no lugar de dados,
### caso tenha problemas com krige.conv()!!!

### Embora o uso da polinomial de 1a ordem tenha representado alguma
melhoria,
### eu não fiquei muito entusiamado por conta do 'lowess' apresentado
### Talvez fosse o caso de testar covariáveis.
### Uma fácil de obter e que pode ter relação é a altitude...
### Mas vamos investigar o efeito nos variogramas

### 3. Distancia máxima e parâmetros para semivariograma
---------------------
hmax <- summary(dados)[[3]][[2]]*.8 ; hmax ### 80% da máxima distância
uvec <- 12; pairs <- 20 ### pairs.min

### 4. Variogramas
-----------------------------------------------------------
### Não precisaria rodar todos, afinal como visto, tá repetindo modelos!!!
### v0:média constante; v1:1st; v2:2nd; v3:cte; v4:coords
v0 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax,
estimator="classical")
v1 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax,
estimator="classical", trend="1st")
v2 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax,
estimator="classical", trend="2nd")
v3 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax,
estimator="classical", trend="cte")
v4 <- variog(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax,
estimator="classical", trend=~coords)

### envelopes
v0.env <- variog.mc.env(dados, obj.var = v0)
v1.env <- variog.mc.env(dados, obj.var = v1)
v2.env <- variog.mc.env(dados, obj.var = v2)
v3.env <- variog.mc.env(dados, obj.var = v3)
v4.env <- variog.mc.env(dados, obj.var = v4)

### 5. Visualizando variogramas
----------------------------------------------
par(par.ori)
par(mfrow=c(2,3))

points(dados, pt.div="quartiles",xlab="Leste", ylab="Norte")

plot(v0, xlab="Dist. (m)", ylim=c(0,summary(v0$v)[6]*1.1), env=v0.env,
main="Isotrópico com tendência")
text(v0$u,v0$v,round(v0$n,1),pos=1);text(v0$u,v0$v,round(v0$u,1),pos=3);lines(v0$u,v0$v);abline(h=mean(v0$v))

plot(v1, xlab="Dist. (m)", ylim=c(0,summary(v1$v)[6]*1.1), env=v1.env,
main="Isotrópico trend 1st")
text(v1$u,v1$v,round(v1$n,1),pos=1);text(v1$u,v1$v,round(v1$u,1),pos=3);lines(v1$u,v1$v);abline(h=mean(v1$v))

plot(v2, xlab="Dist. (m)", ylim=c(0,summary(v2$v)[6]*1.1), env=v2.env,
main="Isotrópico trend 2nd")
text(v2$u,v2$v,round(v2$n,1),pos=1);text(v2$u,v2$v,round(v2$u,1),pos=3);lines(v2$u,v2$v);abline(h=mean(v2$v))

plot(v3, xlab="Dist. (m)", ylim=c(0,summary(v3$v)[6]*1.1), env=v3.env,
main="Isotrópico trend cte")
text(v3$u,v3$v,round(v3$n,1),pos=1);text(v3$u,v3$v,round(v3$u,1),pos=3);lines(v3$u,v3$v);abline(h=mean(v3$v))

plot(v4, xlab="Dist. (m)", ylim=c(0,summary(v4$v)[6]*1.1), env=v4.env,
main="Isotrópico trend coords")
text(v4$u,v4$v,round(v4$n,1),pos=1);text(v4$u,v4$v,round(v4$u,1),pos=3);lines(v4$u,v4$v);abline(h=mean(v4$v))

### Analisando os variogramas
------------------------------------------------
### Realmente o comportamento do histograma de trend='cte' parece conter
uma tendência,
### mas os modelos de 'retirada' usados deixam os dados sem uma estrutura
de
### variação. Dificilmente vai conseguir ajustar um bom modelo nesses
casos, pois
### o comportamento é muito errático. É possível que O erro que você estava
obtendo
### seja decorrente da fragilidade do modelo ajustado.

### Realmente não ficou bom! Mas e aí? O que dá pra fazer? Tenta ssim mesmo?
### Ainda tem que testar a anisotropia, que não tem no seu script original.
?variog4

va <- variog4(dados, uvec=uvec, pairs.min=pairs, max.dist=hmax,
estimator="classical")
par(par.ori)
plot(va)

### Pronto! Apareceu anisotropia...
### Tem duas sugestões de leitura abaixo pra dar uma ajuda...
browseURL('
http://www.rc.unesp.br/igce/aplicada/DIDATICOS/LANDIM/tkrigagem.pdf')
browseURL('
https://www.soils.org/publications/sssaj/abstracts/61/1/SS0610010298')

### Você vai ter que estudar adequadamente a questão da anisotropia em seus
dados
### Isso não te isenta do efeito combinado anisotropia + tendência
### Só pra testar as possibilidades, vamos fazer o variograma na direção
135 graus

### Vai ter que entrar a direções, sendo mais fácil usar PI/n
### Tabela das direções (PI/n) para variogramas
denom <- c(1:4,6,8,10,12,24,36,45,60); k <- (pi/180) ### Denominadores
data.frame(rad=paste0('pi/', denom), rad2=round(pi/denom,4), deg=pi/denom/k)
### 135 é 3*45 ou 3*pi/4

va.135 <- variog(dados, uvec=uvec, pairs.min=pairs, dir=3*pi/4,
max.dist=hmax, estimator="classical")
plot(va.135)
### Aqui parece ser mais fácil ajustar. Talvez tenha que ajustar o 'lag'
### Avaliar a questão da distância entre amostras...

### Lembra do resíduo que separei anteriormente?
### Se fosse usar...
va.r = variog(dados.res, uvec=uvec, pairs.min=pairs, max.dist=hmax,
estimator="classical")
va.r.env <- variog.mc.env(dados.res, obj.var = va.r)
plot(va.r, xlab="Dist. (m)", ylim=c(0,summary(va.r$v)[6]*1.1),
env=va.r.env, main="Residual")

### Vou parar por aqui, porque antes de prosseguir vai precisar ver
### adequadamente essa questão da anisotropia.
### No caso do efeito ser só da anisotropia, vai simplificar a krigagem e
você
### não deve mais ter o problema que vinha tendo.

### <END>

Éder Comunello <c <comunello.eder em gmail.com>omunello.eder em gmail.com>
Dourados, MS - [22 16.5'S, 54 49'W]
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20131026/2b3ca060/attachment.html>