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

Terça Outubro 29 23:25:10 BRST 2013

Caro Eder,

Só hoje vi seu post com CRM anexo.

Muito obrigado pela revisão do CRM.

Como também está a procura de uma maior entendimento da geoestatística,
vamos seguindo.

Vou numerar para separar as observações:

1. Pra entender melhor a área, realmente é café, o recorte esta meio
estranho, mas na realidade são três talhões de 0.5 ha cada.
Cada um tem uma característica bem diferenciada com relação a
produção:alta, baixa  e média.

Segundo li, numa situação em que a média é flutuante, como parece ser o
caso, teríamos de usar krigagem ordinária.

A maioria das observações tinha notado, a linha:
 plot(dados, lowess=T, trend='1st')
achei interessante, mostrando o resultado da retirada de tendencia.

2. no final :

### 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")

Foi possível gerar o semivariograma sem tendencia, com os resíduos, fazer o
ajuste assentimento usando: eyefit(va.r)

Dai fazer a validação, fazendo a validação ou não e gerar o mapa sem
tendencia.

3. Quanto ao envelopamento, não sei se estou certo, mas ele simula valores
máximos e minimos para dependencia, mas como é uma simulação, a cada vez
que vc. roda, pode ter pontos com dependencia, nem que seje um, e as vezes
não dá nenhum, mas de qualquer forma indica que o valor do alcance é
pequeno em relação ao todo e baixa dependencia espacial, como nos meus
dados.

Esse item do envelopamento vale discussão, principalmente em situações que
a variancia sobe além do limite superior e desce, isso significa algo
importante, mas o que?

4. Quanto aos links do Landim já tinha visto, usando o surfer, também vale
atenção.
Fiz o recomendado por ele, mas dobrou os valores do mapa., vc. tentou?

5. Quanto ao ajuste por verosimilhança colocada pelo Prof. Paulo:
Na primeira vez que havia rodado os dados achei muito estranho, pois o
semivariograma subia indefinidamente e o semivariograma teórico por
verossimilhança, usando likfit, baixou muito o patamar.
Interessante se usar likfit com trend="1st", o semivariograma teórico baixa
mais ainda, das dai não consigo krigar

6. Não sei se posso fazer o ajuste assentimento usando o semivariograma com
tendencia mas usar os parametros sem tendencia.

7. Tenho de retirar a tendencia sempre?
Tem casos que há uma patamar claro, mas se rodar o variog com trend, fica
efeito pepita puro...e ai?

8. Outra coisa que não sei se estou entendendo:
após a retirada da tendencia, fazer o mapa com o residuo, mas os valores
não tem nada a ver com os reais, onde estou entendendo errado?

Abraço,

Hélio

Em 26 de outubro de 2013 11:05, Eder Comunello [via R-br] <
ml-node+s2285057n4660724h12 em n4.nabble.com> escreveu:

> 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 <[hidden email]<http://user/SendEmail.jtp?type=node&node=4660724&i=0>[hidden
> email] <http://user/SendEmail.jtp?type=node&node=4660724&i=1>>
> Dourados, MS - [22 16.5'S, 54 49'W]
>
>
>
>
> _______________________________________________
> R-br mailing list
> [hidden email] <http://user/SendEmail.jtp?type=node&node=4660724&i=2>
> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
> código mínimo reproduzível.
>
> ------------------------------
>  If you reply to this email, your message will be added to the discussion
> below:
>
> http://r-br.2285057.n4.nabble.com/Re-R-br-erro-na-validacao-cruzada-com-pacote-geoR-tp4660619p4660724.html
>  To unsubscribe from R-br, click here<http://r-br.2285057.n4.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=3357982&code=aGVsaW9nYWxsb3JvY2hhQGdtYWlsLmNvbXwzMzU3OTgyfC0xMzQ3NTkwMDY4>
> .
> NAML<http://r-br.2285057.n4.nabble.com/template/NamlServlet.jtp?macro=macro_viewer&id=instant_html%21nabble%3Aemail.naml&base=nabble.naml.namespaces.BasicNamespace-nabble.view.web.template.NabbleNamespace-nabble.view.web.template.NodeNamespace&breadcrumbs=notify_subscribers%21nabble%3Aemail.naml-instant_emails%21nabble%3Aemail.naml-send_instant_email%21nabble%3Aemail.naml>
>

-- 
Hélio Gallo Rocha
IFSULDEMINAS - Câmpus Muzambinho
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20131029/7ff7aed5/attachment.html>