rodrigosantana/script_s100.R

## script_s100.R
setwd("k:/2012/geoestatistica/s100")
## Script bastante "automatizado", sem nececissade de incerir dados como: sill
##verificação de dependencia espacial
library(geoR)
library(tcltk)#   para carregar o eyefit, algumas máquinas ele não carrega.
library(scatterplot3d)
## gráficos: primeiro para ver a distribuição(o primeiro gráfico:

## Carregando a base de dados s100
data(s100)

## Atribuindo o objeto geodata s100 à um novo objeto chamado dados
## esta atribuição se deu unica e exclusivamente para não trocar
## o nome do objeto de trabalho ao longo de todo script proposto
## pelo Sr. Hélio originalmente.
dados <- s100

dados##Verificando o objetos criados
n=summary(dados)[[1]] # evita digitar n. de pontos, qaundo necessário
n

##Detecção de dados atípicos
#summary(dados$data)
#var(dados$data)
#cv=((sqrt(var(dados$data)))/(mean(dados$data)))*100
#cv
#hist(dados$data,xlab="Classe de Valores",ylab="Frequência",main="Histograma")
#boxplot(dados$data)
#plot(dados)
#plot(dados, scatter=TRUE, lowes=TRUE) # Gráficos: (dados no espaço), (dados vs coordX), (dados vs coordY), (Z vs Y vs X)

##DETECÇÃO DE TENDÊNCIAS            - como não sabemos a unidade de X1 e X2, convencionei ser metro
## Os dados s100 estão contidos numa área de 1 m x 1 m , acho que nem seria necessário criar uma borda
## mas em casos onde o desenho da área é irregular, teremos de criar uma borda
## para limitar a área da krigagem
limite=read.table("recorte.txt",header=T)
plot(limite,type="l")
points(dados$coords)
points(dados,pt.div="quartile",bor=limite,xlab="Longitude", ylab="Latitude") #   classifica os pontos em quartile

##Distâncias  , normalmente se utiliza 60% da distância máxima
## fiz o calculo para três distâncias, 60%, 80% e 100%, pois cada caso é uma caso

uvec=12     ## pode mudar esse parâmetro
#uvec=seq(0,100,8) ## exemplo de uso do uvec,
##determina pontos de 0 a 100m, com dist.fixa de 8m
pairs.min=30     ## pode mudar esse parâmetro
par=pairs.min
## Plotar quatro gráficas para escolha visual de:
## 1. melhor distância máxima
## 2. Verificar existência de dependencia espacial
par(mfrow=c(1,4))
## hmax com 60 % da dist. máxima
hmax06 = summary(dados)[[3]][[2]]*0.6
##Calculo do semivariograma
vario06= variog(dados,uvec=uvec ,pairs.min=par,max.dist=hmax06, estimator="classical")
#Detecção de dependencia espacial com envelopamento
vario.ent06 <- variog.mc.env(dados, obj.var = vario06)
##Plotando
plot(vario06, envelope = vario.ent06, xl="60% da distância máxima")

## hmax com 80 % da dist. máxima
hmax08 = summary(dados)[[3]][[2]]*0.8
vario08= variog(dados,uvec=uvec,pairs.min=par,max.dist=hmax08, estimator="classical")
vario.ent08 <- variog.mc.env(dados, obj.var = vario08)
plot(vario08, envelope = vario.ent08, xl="80% da distância máxima")

## hmax com 100 % da dist. máxima
hmax01 = summary(dados)[[3]][[2]]
vario01= variog(dados,uvec=uvec,pairs.min=par,max.dist=hmax01, estimator="classical")
vario.ent01 <- variog.mc.env(dados, obj.var = vario01)
plot(vario01, envelope = vario.ent01, xl="100% da distância máxima")

##  escolha da distância para o preparo do semivariograma
hmax=hmax08      ## é esta a distância usada em todo o script
vario= variog(dados,uvec=uvec,pairs.min=par,max.dist=hmax, estimator="classical")
vario.ent <- variog.mc.env(dados, obj.var = vario)
plot(vario, envelope = vario.ent, xl="distância escolhida")

##Significado do variog
##$u = são as distâncias que serão plotadas no eixo no semivariograma.
##$v = é a semivariância para cada distância.
##$n = número de pares de pontos em cada distância.
##$uvec => número de pontos
##sd = desvio padrão naquela distância
##var = variância dos dados
##$n.data = numero de pares trabalhados

##Tabela Resumo das Variancias
#Var06=c(vario06$v)
#Var08=c(vario08$v)
#Var01=c(vario01$v)
#Var=c(vario$v)
#dist06=c(vario06$u)
#dist08=c(vario08$u)
#dist01=c(vario01$u)
#disth=c(vario$u)
#n06=c(vario06$n)
#n08=c(vario08$n)
#n01=c(vario01$n)
#n0=c(vario$n)
#data.frame(n06,dist06,Var06,n08,dist08,Var08,n01,dist01,Var01,n0,disth,Var)
#tabela_resumo=edit(data.frame(n06,dist06,Var06,n08,dist08,Var08,n01,dist01,Var01,n0,disth,Var))

## AJUSTE "À SENTIMENTO"
## Temos de achar um modelo teórico que melhor se adeque ao semivariograma experimental
## obteremos os valores de C0: efeito pepita; alcance: distancia até onde há dep. espacial e
## C1: contribuição
plot(vario, envelope = vario.ent, xl="distância escolhida")
va1=eyefit(vario)
## os valores escolhidos no eyefit que iremos rodar nestas três linhas abaixo,
## serão utilizados nas validações, para escolha do melhor modelo teórico
## achei o modelo gaussiano o mais adequado
sillp=va1[[1]]$cov.pars[1]   ##C1
alcance=va1[[1]]$cov.pars[2] ##phi
pepita=va1[[1]]$nugget    ##C0

va1 ##valores escolhidos , normalmente nem precisa de rodar essas linhas,
    ##é só para verificar se estão certo
sillp
alcance
pepita

##Validação cruzada do modelo "à sentimento"
xv.va1=xvalid(dados,model=va1)##modelo va1

## Pode rodar direto, sem se preocupar, ao final será gerado um relatório
## Este relatório tem a finalidade de ajudar na escolha do melhor modelo
## alguns dados podem ser retirados do relatório
mstdva1=sum(xv.va1$std.error)/n##erro médio reduzido
sqrtva1=sqrt((sum(xv.va1$std.error^2))/n)##Desvio Padrão dos erros reduzidos
GDEva1= (sillp/(sillp+pepita))*100
corva1=cor(xv.va1$data,xv.va1$predicted)
reglin_va1=lm(xv.va1$data~xv.va1$predicted)
##data.frame(xv.va1$data,xv.va1$predicted, xv.va1$error, xv.va1$std.error, xv.va1$krige.var)
##edit(data.frame(xv.va1$data,xv.va1$predicted, xv.va1$error, xv.va1$std.error, xv.va1$krige.var))
##plot(reglin_va1)
summary(reglin_va1)
r2_va1=summary(reglin_va1)[[8]]
r2_va1
merrova1=mean(xv.va1$error)##média do erro de VA1 para compor a tebela
reglin_va1  ## valores da equação de regressão
summary(reglin_va1)
plot(xv.va1$data,xv.va1$predicted)
abline(reglin_va1)
summary(reglin_va1)[[8]]   ## s.squared
anova(reglin_va1)

##MÉTODO DA VEROSSIMILHANÇA >>>> ML
ml=likfit(dados,ini=c(sillp,alcance),nugget=pepita)##C(SIGMASQ,PHI)##MÉTODO DA VEROSSIMILHANÇA>ML???
#ml
tausqml=ml$tausq ##C0 pepipa
sigmasqml= ml$sigmasq	##C1 patamar
phiml=ml$phi ## alcance
plot(vario,xlab="Distancia (h)",ylab="Semivariancia")
lines(ml,col="red")
xv.ml<-xvalid(dados,model=ml)
##xv.ml
mstdml=sum(xv.ml$std.error)/n   ## média quanto mais próximo de 0 melhor
sqrtml=sqrt((sum(xv.ml$std.error^2))/n)  ## desvio padrão quanto mais próximo de 1
GDEml = (sigmasqml/(sigmasqml+tausqml))*100
corml=cor(xv.ml$data,xv.ml$predicted)
reglin_ml=lm(xv.ml$data~xv.ml$predicted)
r2_ml=summary(reglin_ml)[[8]]
r2_ml
merroml=mean(xv.ml$error)
reglin_ml
anova(reglin_ml)
plot(xv.ml$data,xv.ml$predicted,xlab="Dados",ylab="Prediicted_ML", main="Gráfico de Validação Cruzada Verossimilhança \n  Dados X Predição",cex.main=1.3)
abline(reglin_ml)
summary(reglin_ml)[[8]]  ##r2

##QUADRADOS MINIMOS ORDINAL >>>>0LS
##modelo de validação cruzada
ols=variofit(vario,ini=c(sillp,alcance),nugget=pepita,wei="equal")
plot(vario,xlab="Distancia (h)",ylab="Semivariancia")
lines(ols,col="red")
xv.ols=xvalid(dados,model=ols)
##xv.ols
#ols
tauols=ols[[1]]   	##C0
sigmaols =ols[[2]][1]	##C1
phiols=ols[[2]][2]	##alcance
mstdols=sum(xv.ols$std.error)/n
sqrtols=sqrt((sum(xv.ols$std.error^2))/n)
GDEols = (sigmaols/(sigmaols+tauols))*100
reglin_ols=lm(xv.ols$data~xv.ols$predicted)
r2_ols=summary(reglin_ols)[[8]]
r2_ols
merrools=mean(xv.ols$error)
reglin_ols
anova(reglin_ols)
plot(xv.ols$data,xv.ols$predicted)
abline(reglin_ols)
corols=cor(xv.ols$data,xv.ols$predicted)
summary(reglin_ols)[[8]]  ##r2
##QUADRADOS MINIMOS PONDERADO >>> WLS????
wls=variofit(vario,ini=c(sillp,alcance),nugget=pepita)
plot(vario,xlab="Distancia (h)",ylab="Semivariancia", main="Quadrados Minímos Ponderados")
lines(wls,col="red")
xv.wls=xvalid(dados,model=wls)
#xv.wls
#wls
tauwls=wls[[1]]   	##C0
sigmawls =wls[[2]][1]	##C1
phiwls=wls[[2]][2]	##alcance
mstdwls=sum(xv.wls$std.error)/n
sqrtwls=sqrt((sum(xv.wls$std.error^2))/n)
GDEwls = (sigmawls/(sigmawls+tauwls))*100
corwls=cor(xv.wls$data,xv.wls$predicted)
reglin_wls=lm(xv.wls$data~xv.wls$predicted)
r2_wls=summary(reglin_wls)[[8]]
r2_wls
merrowls=mean(xv.wls$error)
reglin_wls
anova(reglin_wls)
plot(xv.wls$data,xv.wls$predicted)
abline(reglin_wls)
summary(reglin_wls)#[[8]]  ##r2
##MODELO GAUSSIANO
wlsg=variofit(vario,ini=c(sillp,alcance),nugget=pepita,cov.model="gaus")##C(SIGMASQ,PHI)
plot(vario,xlab="Distancia (h)",ylab="Semivariancia", main="Modelo Gaussiano")
lines(wlsg,col="maroon")
xv.wlsg=xvalid(dados,model=wlsg)
##xv.wlsg
#wlsg
tauwlsg=wlsg[[1]]   	##C0
sigmawlsg =wlsg[[2]][1]	##C1
phiwlsg=wlsg[[2]][2]	##alcance
mstdwlsg=sum(xv.wlsg$std.error)/n
sqrtwlsg=sqrt((sum(xv.wlsg$std.error^2))/n)
GDEwlsg = (sigmawlsg/(sigmawlsg+tauwlsg))*100
corwlsg=cor(xv.wlsg$data,xv.wlsg$predicted)
corwlsg
reglin_wlsg=lm(xv.wlsg$data~xv.wlsg$predicted)
reglin_glm_wlsg=glm(xv.wlsg$data~xv.wlsg$predicted)
reglin_glm_wlsg
reglin_wlsg
r2_wlsg=summary(reglin_wlsg)[[8]] ##fator de deteminação da validação curuzada entre data X predicted
r2_wlsg *100
merrowlsg=mean(xv.wlsg$error)
reglin_wlsg
anova(reglin_wlsg)
plot(xv.wlsg$data,xv.wlsg$predicted)
abline(reglin_wlsg)
summary(reglin_wls)[[8]]  ##r2


##Cálculo do Erro Absoluto    --- Somatória???
ea_va1=sum((xv.va1$predicted)-(xv.va1$data))
ea_ml=sum((xv.ml$predicted)-(xv.ml$data))
ea_ols=sum((xv.ols$predicted)-(xv.ols$data))
ea_wls=sum((xv.wls$predicted)-(xv.wls$data))
ea_wlsg=sum((xv.wlsg$predicted)-(xv.wlsg$data))

##significado da validação cruzada xvalid
#?xvalid
#$data - os dados originais.
#$predicted - os valores previstos pela validação cruzada.
#$krige.var	- a variação de previsão de validação cruzada.
#$error	- as diferenças de dados - valor previsto ($data-$predicted) .
#$std.error  os erros dividida pela raiz quadrada das variâncias de predição.
#$prob - 	a probabilidade cumulativa no valor original com uma distribuição normal com parâmetros dados com os resultados de validação cruzada.
#?variog
##Significado do variog   vario
##$u = são as distâncias que serão plotadas no eixo no semivariograma.
##$v = é a semivariância para cada distância.
##$n = número de pares de pontos em cada distância.
##$uvec => número de pontos
##sd = desvio padrão naquela distância
##var = variância dos dados
##$n.data = numero de pares trabalhados.

##PLOTANDO TODOS OS MODELOS
plot(vario,xlab="Distancia (h)",ylab="Semivariancia")
lines(va1,col="red")
lines(ml,col="green")
lines(ols,col="maroon")
lines(wls, col="black")
lines(wlsg,col="cyan")


 summary(dados)
 ###KRIGAGEM
 ## serão calculadas cinco krigagens correspondendo aos cinco modelos de ajuste
loci=expand.grid(seq(-0.0391304,1.0391304,l=400),seq(-0.0391304,1.0391304,l=400))
limite=read.table("recorte.txt",header=T)

kc=krige.conv(dados,loc=loci,border=limite,krige=krige.control(obj=va1))
kc2=krige.conv(dados,loc=loci,border=limite,krige=krige.control(obj=ml))
kc3=krige.conv(dados,loc=loci,border=limite,krige=krige.control(obj=ols))
kc4=krige.conv(dados,loc=loci,border=limite,krige=krige.control(obj=wls))
kc5=krige.conv(dados,loc=loci,border=limite,krige=krige.control(obj=wlsg))

l=4#escolha do n. de classificações do mapa, pode mudar

##Todas as Krigagens, não consegui ajustar exatamente os pontos com a krigagem, ficou um pouco deslocado.
par(mfrow=c(2,3))
points(dados,pt.div="quartile",bor=limite,  xlab="Longitude", ylab="Latitude")
image(kc,loc=loci,border=limite,val=kc$predict,col="gray"(seq(1,0.1,l=l)),xlab=" ",ylab=" ", main="Krigagem à sentimento")
#legend.krige(c(0.2,0.8),c(1,1.1),cex.leg=0.7,val=kc2$predict,col=gray(seq(1,0.1,l=l)),off=0.8)
par(new=TRUE)
points(dados,pt.div="quartile",bor=limite,  xlab="Longitude", ylab="Latitude")
#points(dados$coords)
image(kc2,loc=loci,border=limite,val=kc2$predict,col="gray"(seq(1,0.1,l=l)),xlab=" ",ylab=" ", main="Krigagem com Modelo Verosimilhança - ML")
#legend.krige(c(0.2,0.8),c(-0.2,-0.1),cex.leg=0.7,val=kc2$predict,col=gray(seq(1,0.1,l=l)),off=0.8)
par(new=TRUE)
points(dados,pt.div="quartile",bor=limite)#,xlab="Longitude", ylab="Latitude")
#points(dados$coords)
image(kc3,loc=loci,border=limite,val=kc3$predict,col="gray"(seq(1,0.1,l=l)),xlab=" ",ylab=" ", main="Krigagem com Modelo Mínimo Ordinal- OLS")
#legend.krige(c(0.2,0.8),c(-0.2,-0.1),cex.leg=0.7,val=kc3$predict,col=gray(seq(1,0.1,l=l)),off=0.8)
par(new=TRUE)
points(dados,pt.div="quartile",bor=limite,xlab="Longitude", ylab="Latitude")
#points(dados$coords)
image(kc4,loc=loci,border=limite,val=kc4$predict,col="gray"(seq(1,0.1,l=l)),xlab=" ",ylab=" ", main="Krigagem com Modelo Mínimo Ponderado- WLS")
#legend.krige(c(0.2,0.8),c(-0.2,-0.1),cex.leg=0.7,val=kc4$predict,col=gray(seq(1,0.1,l=l)),off=0.8)
par(new=TRUE)
points(dados,pt.div="quartile",bor=limite,xlab="Longitude", ylab="Latitude")
#points(dados$coords)
image(kc5,loc=loci,border=limite,val=kc5$predict,col="gray"(seq(1,0.1,l=l)),xlab=" ",ylab=" ", main="Krigagem com Modelo Gaussiano- WLSG")
#legend.krige(c(0.2,0.8),c(-0.2,-0.1),cex.leg=0.7,val=kc4$predict,col=gray(seq(1,0.1,l=l)),off=0.8,)
par(new=TRUE)
points(dados,pt.div="quartile",bor=limite,xlab="Longitude", ylab="Latitude")

## perceberam que tentei colocar no mesmo mapa o mapa da krigagem e os pontos classificados em quartile
## ficou meio deslocado, quem puder ajudar a ajustar...
## tenho um método, que é bastante chato, quem souber com autoajustar, passe...
## se ninguém conseguir, eu passo como faço, ai está o começo
## 1. plotar primeiro os pontos classificados em quartile
 #points(dados,pt.div="quartile",bor=limite)
 #2. Meço no monitor a distãncia da borda do ponto (0,0), no meu deu 0,5 cm e
 ## de (0, 0,2) deram 2,25 cm  , faço uma regra de três..........


 ######    Montagem da tabela comparativa de resultados
##Nome dos Modelos
 Asentimento=c("Asentimento")
 Verossimilhança=c("Verossimilhança-ML")
 OLS=c("Minimo Ordinal-OLS")
 WLS=c("Minimo Ponderado-WLS")
 WLSG=c("Gaussiano")
 Minimo=c("Minimo")
 Máximo=c("Máximo")
 MODELO=c(Asentimento,Verossimilhança,OLS,WLS,WLSG,Minimo,Máximo)


 c0=c(pepita,tausqml,tauols,tauwls,tauwlsg)
 c0=c(pepita,tausqml,tauols,tauwls,tauwlsg,summary(c0)[[1]],summary(c0)[[6]])
 c1=c(sillp,sigmasqml,sigmaols,sigmawls,sigmawlsg)
 c1=c(sillp,sigmasqml,sigmaols,sigmawls,sigmawlsg,summary(c1)[[1]],summary(c1)[[6]])
 Patamar=c(c0+c1)
 Alcance=c(alcance,phiml,phiols,phiwls,phiwlsg)
 Alcance=c(alcance,phiml,phiols,phiwls,phiwlsg,summary(Alcance)[[1]],summary(Alcance)[[6]])
 ER=c(mstdva1,mstdml,mstdols,mstdwls,mstdwlsg)    ##erro médio reduzido
 ER=c(mstdva1,mstdml,mstdols,mstdwls,mstdwlsg,summary(ER)[[1]],summary(ER)[[6]])##erro médio reduzido
 mean_error=c(merrova1,merroml,merrools,merrowls,merrowlsg)#Média do Erro
 mean_error=c(merrova1,merroml,merrools,merrowls,merrowlsg,summary(mean_error)[[1]],summary(mean_error)[[6]])#Média do Erro
 Ser=c(sqrtva1,sqrtml,sqrtols,sqrtwls,sqrtwlsg)##Desvio Padrão dos erros reduzidos
 Ser=c(sqrtva1,sqrtml,sqrtols,sqrtwls,sqrtwlsg,summary(Ser)[[1]],summary(Ser)[[6]])##Desvio Padrão dos erros reduzidos
 SerM=c(2-sqrtva1,2-sqrtml,2-sqrtols,2-sqrtwls,2-sqrtwlsg)
 SerM=c(2-sqrtva1,2-sqrtml,2-sqrtols,2-sqrtwls,2-sqrtwlsg,2-summary(SerM)[[1]],2-summary(SerM)[[6]])
 GDE=c(GDEva1,GDEml,GDEols,GDEwls,GDEwlsg)
 GDE=c(GDEva1,GDEml,GDEols,GDEwls,GDEwlsg,summary(GDE)[[1]],summary(GDE)[[6]])
 R2_dataXpredictede=c(r2_va1,r2_ml,r2_ols,r2_wls,r2_wlsg)
 R2_dataXpredictede=c(r2_va1,r2_ml,r2_ols,r2_wls,r2_wlsg,summary(R2_dataXpredictede)[[1]],summary(R2_dataXpredictede)[[6]])
 Correlação=c(corva1,corml,corols,corwls,corwlsg)
 Correlação=c(corva1,corml,corols,corwls,corwlsg,summary(Correlação)[[1]],summary(Correlação)[[6]])
 EA=c(ea_va1,ea_ml,ea_ols,ea_wls,ea_wlsg)
 EA=c(ea_va1,ea_ml,ea_ols,ea_wls,ea_wlsg,summary(EA)[[1]],summary(EA)[[6]])
mean_krigevar=c(summary(kc$krige.var)[[4]],summary(kc2$krige.var)[[4]],summary(kc3$krige.var)[[4]],summary(kc4$krige.var)[[4]],summary(kc5$krige.var)[[4]])
mean_krigevar=c(summary(kc$krige.var)[[4]],summary(kc2$krige.var)[[4]],summary(kc3$krige.var)[[4]],summary(kc4$krige.var)[[4]],summary(kc5$krige.var)[[4]],summary(mean_krigevar)[[1]],summary(mean_krigevar)[[6]])
mean_predict=c(summary(kc$predict)[[4]],summary(kc2$predict)[[4]],summary(kc3$predict)[[4]],summary(kc4$predict)[[4]],summary(kc5$predict)[[4]])
mean_predict=c(summary(kc$predict)[[4]],summary(kc2$predict)[[4]],summary(kc3$predict)[[4]],summary(kc4$predict)[[4]],summary(kc5$predict)[[4]],summary(mean_predict)[[1]],summary(mean_predict)[[6]])


## veja os valores AIC, escolha o menor e o maior
AIC(reglin_va1) ##reglin=regressão linear ...reglin_av1=lm(xv.va1$data~xv.va1$predicted)   , sendo va1 - asentimento
AIC(reglin_ml)
AIC(reglin_ols)
AIC(reglin_wls)
AIC(reglin_wlsg)
##especificar aqui o AIC menor e o maior
aic_min=(AIC(reglin_ml))##especificar
aic_max=(AIC(reglin_wlsg))##especificar
aic_min
aic_max
dif_va1=((AIC(reglin_va1))-(aic_min))
dif_ml=((AIC(reglin_ml))-(aic_min))
dif_ols=((AIC(reglin_ols))-(aic_min))
dif_wls=((AIC(reglin_wls))-(aic_min))
dif_wlsg=((AIC(reglin_wlsg))-(aic_min))
dif_va1
dif_ml
dif_ols
dif_wls
dif_wlsg
AKAIKE=c(dif_va1,dif_ml,dif_ols,dif_wls,dif_wlsg,aic_min,aic_max)
AK=data.frame(AKAIKE)
#BIC(reglin_va1)
#BIC(reglin_ml)
#BIC(reglin_ols)
#BIC(reglin_wls)
#BIC(reglin_wlsg)
#bic_min=(BIC(reglin_ml))##especificar
#bic_max=(BIC(reglin_wlsg))##especificar
##bic_min
#bic_max
#difb_va1=((BIC(reglin_va1))-(bic_min))
#difb_ml=((BIC(reglin_ml))-(bic_min))
#difb_ols=((BIC(reglin_ols))-(bic_min))
#difb_wls=((BIC(reglin_wls))-(bic_min))
##difb_wlsg=((BIC(reglin_wlsg))-(bic_min))
##difb_va1
##difb_ml
##difb_ols
##difb_wls
##difb_wlsg
##BIC=c(difb_va1,difb_ml,difb_ols,difb_wls,difb_wlsg,bic_min,bic_max)
##data.frame(BIC)
tabela_resumo_B=edit(data.frame(MODELO,c0,c1,Patamar,Alcance, ER, Ser,SerM,mean_error,GDE,R2_dataXpredictede,Correlação,EA,mean_krigevar,mean_predict,AKAIKE))#,BIC))
write.table(tabela_resumo_B, file="k:\\2012/geoestatistica/s100/akaike.xls", sep = '\t',col.names = NA)


# Variogramas direcionais:
vario.4 <- variog4(dados, max.dist = hmax)
plot(vario.4, lwd=3)
plot(variog(dados, max.dist=hmax))
lines.variomodel(cov.model="exp", cov.pars=c(1,.3), nug=0, max.dist=hmax)
lines.variomodel(cov.model="mat", cov.pars=c(.85,.2), nug=0.1, kappa=1,max.dist=hmax, lty=2)
lines.variomodel(cov.model="sph", cov.pars=c(.8,.8), nug=0.1,max.dist=hmax, lwd=2)


par(mfrow=c(2,3))
vario0=variog(dados,max.dist=hmax,pairs.min=par,direction=pi,estimator="modulus")
vario45 = variog(dados,max.dist=hmax, pairs.min=par,direction=pi/4,estimator="modulus")
vario90 = variog(dados,max.dist=hmax, pairs.min=par,direction=pi/2,estimator="modulus")
vario135 = variog(dados,max.dist=hmax, pairs.min=par,direction=   (pi/4 + pi/2) ,estimator="modulus")
plot(vario0, xlab="zero")
plot(vario45, xlab="quarenta e cinco")
plot(vario90, xlab="noventa")
plot(vario180, xlab="cento e trinta e cinco")

###########################
#Após escolha da direção, fazer ajuste, no caso foi ? graus
variodir=vario45
variodir = variog(dados,max.dist=hmax,direction=pi/4,estimator="modulus")
plot(variodir)
dir=eyefit(variodir) # Ajusta a olho (a sentimento) (chute inicial para o método iterativo (método numérico))
summary (dir)
taudir=dir[[1]]$nugget 	 ##C0
sigmasqdir =dir[[1]]$cov.pars[1]	##C1
phidir=dir[[1]]$cov.pars[2] ##phi
#####valores escolhidos
taudir
sigmasqdir
phidir

xv.dir=xvalid(dados,model=dir)##modelo modulous com direção escolhida
xv.dir
mean(xv.dir$error)
mean(xv.dir$std.error)
sqrt((sum(xv.dir$std.error^2))/n)
GDEdir=(sigmasqdir/(sigmasqdir+tausqdir))*100
GDEdir

plot(xv.dir$data, xv.dir$predict)
sqrt((sum(xv.dir$predict^2))/n)
sqrt((sum(xv.dir$error^2))/n)

wlsdir=variofit(variodir,ini=c(sigmasqdir,phidir),nugget=tausqdir)
plot(variodir,xlab="Distancia (h)",ylab="Semivariancia")
lines(wls,col="red")
xv.wlsdir=xvalid(dados,model=wls)
xv.wlsdir
wlsdir
tauwlsdir=wlsdir[[1]]   	##C0
sigmawlsdir =wlsdir[[2]][1]	##C1
phiwlsdir=wlsdir[[2]][2]	##alcance
#####valores escolhidos
tauwlsdir
sigmawlsdir
phiwlsdir

plot(xv.wlsdir$data,xv.wlsdir$predicted)
cor(xv.wlsdir$data,xv.wlsdir$predicted)

sum(xv.wlsdir$std.error)/n
sqrt((sum(xv.wlsdir$std.error^2))/n)
GDEwls = (sigmawlsdir/(sigmawlsdir+tauwlsdir))*100
GDEwls
plot(variodir, xlab="Distância escolhida(h)", ylab="Semivariograma (classical)")
lines(wlsdir)


##$u = são as distâncias que serão plotadas no eixo no semivariograma.
##$v = é a semivariância para cada distância.
##$n = número de pares de pontos em cada distância.
##$uvec => número de pontos
##sd = desvio padrão naquela distância
##var = variância dos dados
##$n.data = numero de pares trabalhados.


####vario2 - modulus
vario2 = variog(dados, max.dist=hmax, pairs.min=30, estimator="modulus") # estimador do semivariograma (estimador robusto de Cressie e Hawkins, 1980), faz o semivariograma considerando o estimador robusto, dados sem outliers os dois ficam parecidos, qdo fica diferente eles ficam muito diferentes...
vario2
vario2$ind.bin
plot(vario2, xlab="Distância (h)", ylab="Semivariograma (Cressie e Hawkins)")#SÓ PARA DADOS COM OUTLIERS
va2=eyefit(vario2)##esférico???
va2#####valores escolhidos
sigmava2=va2[[1]]$cov.pars[1]   ##C1
phiva2=va2[[1]]$cov.pars[2] ##phi
tauva2=va2[[1]]$nugget 	 ##C0

sigmava2
phiva2
tauva2

plot(xv.va2$data,xv.va2$predicted)
cor(xv.va2$data,xv.va2$predicted)

xv.va2=xvalid(dados,model=va2)##modelo va2 = sph?
xv.va2   ## modelo exponencial
sum(xv.va2$std.error)/n
sqrt((sum(xv.va2$std.error^2))/n)
GDEva2= (sigmava2/(sigmava2+tauva2))*100
GDEva2
gde1=(tauva2/(sigmava2+tauva2))*100
gde1

plot(vario2, xlab="Distância (h)", ylab="Semivariograma (Cressie e Hawkins)")
lines(va2,col="red")

summary(xv.va2)
summary(dados)


anova((data.frame(MODELO,c0,c1,Patamar,Alcance,mean_std_error,mean_error,SQRT,SQRTm,GDE,R2_dataXpredictede,Correlação)))##data.frame(nomes,mean_krigevar,mean_predict))


summary(kc$krige.var)[[4]]
summary(kc$predict)##[[4]]

###KRIGAGEM VA1
boxplot(kc$krige.var)
boxplot(kc$predict)
hist(kc$krige.var)###data,xlab="Classe de Valores",ylab="Frequência",main="Histograma")
hist(kc$predict)


###KRIGAGEM ml
kc2=krige.conv(dados,loc=loci,border=limite,krige=krige.control(obj=ml))
image(kc2,loc=loci,border=limite,val=kc2$predict,col="gray"(seq(1,0.1,l=5)))##l=5(n. classes)
legend.krige(c(348030,348210),c(7637910,7637925),cex.leg=0.7,val=kc2$predict,col=gray(seq(1,0.1,l=5)),off=0.8)


summary(kc2$krige.var)[[4]]
summary(kc2$predict)[[4]]

boxplot(kc2$krige.var)
boxplot(kc2$predict)
hist(kc2$krige.var)
hist(kc2$predict)

###KRIGAGEM ols
kc3=krige.conv(dados,loc=loci,border=limite,krige=krige.control(obj=ols))
image(kc3,loc=loci,border=limite,val=kc3$predict,col="gray"(seq(1,0.1,l=5)))##l=5(n. classes)
legend.krige(c(348030,348210),c(7637910,7637925),cex.leg=0.7,val=kc3$predict,col=gray(seq(1,0.1,l=5)),off=0.8)

summary(kc3$krige.var)[[4]]
summary(kc3$predict)[[4]]


boxplot(kc3$krige.var)
boxplot(kc3$predict)
hist(kc3$krige.var)###data,xlab="Classe de Valores",ylab="Frequência",main="Histograma")
hist(kc3$predict)

###KRIGAGEM wls
kc4=krige.conv(dados,loc=loci,border=limite,krige=krige.control(obj=wls))
image(kc4,loc=loci,border=limite,val=kc4$predict,col="gray"(seq(1,0.1,l=5)))##l=5(n. classes)
legend.krige(c(348030,348210),c(7637910,7637925),cex.leg=0.7,val=kc4$predict,col=gray(seq(1,0.1,l=5)),off=0.8)

summary(kc4$krige.var)[[4]]
summary(kc4$predict)[[4]]

boxplot(kc4$krige.var)
boxplot(kc4$predict)
hist(kc4$krige.var)###data,xlab="Classe de Valores",ylab="Frequência",main="Histograma")
hist(kc4$predict)

###KRIGAGEM wlsg
kc5=krige.conv(dados,loc=loci,border=limite,krige=krige.control(obj=wlsg))
image(kc5,loc=loci,border=limite,val=kc5$predict,col="gray"(seq(1,0.1,l=5)))##l=5(n. classes)
legend.krige(c(348030,348210),c(7637910,7637925),cex.leg=0.7,val=kc4$predict,col=gray(seq(1,0.1,l=5)),off=0.8)

summary(kc4$krige.var)[[4]]
summary(kc4$predict)[[4]]

boxplot(kc4$krige.var)
boxplot(kc4$predict)
hist(kc4$krige.var)###data,xlab="Classe de Valores",ylab="Frequência",main="Histograma")
hist(kc4$predict)