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] < [email protected]> 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 [email protected] 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.
