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.,Em 29 de outubro de 2013 22:25, Hélio Gallo Rocha <heliogallorocha@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@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)### 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 geralpar(mfrow=c(1,1)); par()$mfrowpoints(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.geodataplot(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ômaloplot(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.spatialhead(unclass(trend.spatial('cte', dados))) ### constantehead(unclass(trend.spatial(~coords, dados))) ### polinomal 1a ordem, usando coordshead(unclass(trend.spatial('1st', dados))) ### mesmo que anteriorhead(unclass(trend.spatial('2nd', dados))) ### polinomal 2a ordem, usando coordshead(unclass(trend.spatial(~coords[,1], dados))) ### em função de Xhead(unclass(trend.spatial(~coords[,2], dados))) ### em função de Yhead(unclass(trend.spatial(~coords+I(coords[,1]^2), dados))) ### coords + X**2head(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 2ndhead(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)$residualstrend.fit <- lm(dados$data ~ trend.matrix)$fitted.values### 'Empacotando' os resíduos em um geodata-----------------------------------dados.res <- dadosnames(trend.res) <- NULL; trend.resdados.res$data <- trend.res; dados.res$datastr(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ânciauvec <- 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:coordsv0 <- 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)### envelopesv0.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.?variog4va <- 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...### 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 variogramasdenom <- c(1:4,6,8,10,12,24,36,45,60); k <- (pi/180) ### Denominadoresdata.frame(rad=paste0('pi/', denom), rad2=round(pi/denom,4), deg=pi/denom/k)### 135 é 3*45 ou 3*pi/4va.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>
_______________________________________________
R-br mailing list
[hidden email]
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.
http://r-br.2285057.n4.nabble.com/Re-R-br-erro-na-validacao-cruzada-com-pacote-geoR-tp4660619p4660724.htmlIf you reply to this email, your message will be added to the discussion below:--
Hélio Gallo Rocha
IFSULDEMINAS - Câmpus Muzambinho
_______________________________________________
R-br mailing list
R-br@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.