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

Quinta Outubro 31 16:21:03 BRST 2013

Caro Eder,

Já havia feito a retirada de tendencia com o atributo altura da planta, que
está ligado com a produção.

O variograma gerado ficou diferente se usar polinomial de 1st, mas não sei
como validar isso depois...

Abraço

Em 31 de outubro de 2013 13:15, Éder Comunello
<comunello.eder em gmail.com>escreveu:

> Boa tarde, Helio!
>
> Vi seu email, mas vou precisar reservar um tempo pra tratar dele. Você
> chegou a usar as glebas como covariável na modelagem? Para isso cada ponto
> deveria ter uma var a mais identificando a gleba origem.
>
> Atte.,
>
>
>
> Éder Comunello <c <comunello.eder em gmail.com>omunello.eder em gmail.com>
> Dourados, MS - [22 16.5'S, 54 49'W]
>
>
> Em 29 de outubro de 2013 22:25, Hélio Gallo Rocha <
> heliogallorocha em gmail.com> escreveu:
>
>> 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
>>
>> _______________________________________________
>> R-br mailing list
>> R-br em listas.c3sl.ufpr.br
>> 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.
>>
>
>

-- 
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/20131031/d862a3b9/attachment.html>