Apresentação Espacial

Introdução

O banco de dados deste trabalho está disponível no software Rstudio, “head” que traz em seu interior 29 observações de aquiferos, contendo coodenadas “x, y” e os dados em questão “z” em determinada região de estudo.

Vamos aos cálculos

Para da início as atividades, se necessário fazer a istalação dos pacotes (geoR) e (MASS), caso não se tenha previamente instalado em sua maquina, e a requerição dos mesmos.

if(!require(geoR)){install.packages("geoR")}

## Loading required package: geoR

## Warning: package 'geoR' was built under R version 3.4.4

## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.7-5.2.1 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------

if(!require(MASS)){install.packages("MASS")}

## Loading required package: MASS

## Warning: package 'MASS' was built under R version 3.4.4

Sobre os Dados

Os dados disponíveis no software Rstudio foram reorganizados e salvos em formato “txt”, afim de se ter uma melhor visualização das informações.

setwd("F:/Estatistica Espacial")
dados <- read.table("BDinf.txt", head=T); dados

##        x     y    z
## 1   6.86  6.41 1061
## 2   4.32  5.02 1194
## 3   5.98  6.01 1117
## 4  11.61  4.99  880
## 5   5.52  3.79 1202
## 6  10.87  8.27  757
## 7   8.61  3.92 1038
## 8  12.64  6.77  817
## 9  14.70 10.43  630
## 10 13.91 10.91  617
## 11  9.47  5.62  986
## 12 14.36 11.03  625
## 13  8.99  7.31  840
## 14 11.93  6.78  847
## 15 11.75 10.80  645
## 16 13.23 10.18  662
## 17 13.56  9.74  685
## 18  8.06  5.76 1023
## 19 10.95  3.72  998
## 20 14.71 11.41  584
## 21 16.27  7.27  611
## 22 12.33  6.87  847
## 23 13.01  7.05  745
## 24 13.56  7.42  725
## 25 13.76  8.35  688
## 26 12.54  9.04  676
## 27  8.97  8.60  768
## 28  9.22  8.55  782
## 29  9.64  3.38 1022

Informações dos dados

Leitura dos dados

head(dados)

##       x    y    z
## 1  6.86 6.41 1061
## 2  4.32 5.02 1194
## 3  5.98 6.01 1117
## 4 11.61 4.99  880
## 5  5.52 3.79 1202
## 6 10.87 8.27  757

Dimensão dos dados

dim(dados)

## [1] 29  3

Nome das colunas dos dados

colnames(dados)

## [1] "x" "y" "z"

Quantidade de observações dos dados

rownames(dados)

##  [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14"
## [15] "15" "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28"
## [29] "29"

Dimenção dobanco de dados

dimnames(dados)

## [[1]]
##  [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12" "13" "14"
## [15] "15" "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28"
## [29] "29"
## 
## [[2]]
## [1] "x" "y" "z"

Análise exploratoria dos dados

Será dado o nome de “INF” para informação dos dados.

INF <-as.geodata(dados)

Resumo dos dados

Através do comando Summary, temos uma série de informações sobre os dados, Ou seja, o valor maximo e minimo das coodenadas, a distância minima e máxima entre as observações, e média e etc…

summary(INF)

## Number of data points: 29 
## 
## Coordinates summary
##         x     y
## min  4.32  3.38
## max 16.27 11.41
## 
## Distance summary
##       min       max 
##  0.254951 12.197713 
## 
## Data summary
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  584.000  676.000  782.000  830.069  998.000 1202.000

visualização gráfica dos dados

Gráfico dos dados originais sem adição de efeito, para visualizar seu respectivo comportamento.

plot(INF)

Gráfico dos dados originais sem adição de efeito, com a linha de ajuste de estacionáriedade.

plot(INF, low=T)

Gráfico dos dados com efeito de tendência de primeira ordem “1st”

plot(INF, low=T, trend="1st")

Adição de algumas perfumarias para melhorar a estética do gráfico dos dados com efeito de tendência de primeira ordem “1st”, com adição dos quintis.

args(points.geodata)

## function (x, coords = x$coords, data = x$data, data.col = 1, 
##     borders = x$borders, pt.divide = c("data.proportional", "rank.proportional", 
##         "quintiles", "quartiles", "deciles", "equal"), lambda = 1, 
##     trend = "cte", abs.residuals = FALSE, weights.divide = "units.m", 
##     cex.min, cex.max, cex.var, pch.seq, col.seq, add.to.plot = FALSE, 
##     x.leg, y.leg = NULL, dig.leg = 2, round.quantiles = FALSE, 
##     permute = FALSE, ...) 
## NULL

points(INF, trend="1st", pt.div="quint", perm=T)

Variograma

Variograma dos dados originais sem adição de efeito e seu respectivo gráfico.

INFO <- variog(INF)

## variog: computing omnidirectional variogram

plot(INFO)

Variograma com efeito de tendência de primeira ordem “1st” e seu respectivo gráfico.

INFO <- variog(INF, trend = "1st")

## variog: computing omnidirectional variogram

plot(INFO)

Variograma com efeito de tendência na coordenada X e seu respectivo gráfico.

INFO <- variog(INF,trend = ~coords[,1])

## variog: computing omnidirectional variogram

plot(INFO)

Variograma com efeito de tendência na coordenada y e seu respectivo gráfico.

INFO <- variog(INF,trend = ~coords[,2])

## variog: computing omnidirectional variogram

plot(INFO)

Variograma com 1/3 da distância máxima, e seu respectivo gráfico.

INFO <- variog(INF, max.dist=5)

## variog: computing omnidirectional variogram

plot(INFO)

O variograma com efeito de tendência de primeira ordem “1st”, mostra aparentemente um comportamento com maior influência com efeito de tendencia de primeira ordem.

INFO1 <- variog(INF, trend = "1st")

## variog: computing omnidirectional variogram

plot(INFO1)

Envelope simulado

Gráfico do envelope simulado para verificar se existe pontos fora dos limites de predição “Envelope”.

INFO1.env <- variog.mc.env(INF, obj=INFO)

## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops

plot(INFO, env=INFO1.env)
lines(INFO1.env)

Máxima Verosimilhança

Ajuste dos modelos pelo método de máxima verosimilhança, afim de encontrar um melhor modelo que se adeque com os dados.

Função wave

Ajuste do modelo inicial “ml1” sem efeito de tendência e sem trasformação.

ml1 <- likfit(INF, ini=c(343145.5, 6.66 ),nug= 0)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml1

## likfit: estimated model parameters:
##        beta       tausq     sigmasq         phi 
## "   865.85" "     0.00" "103420.24" "    79.81" 
## Practical Range with cor=0.05 for asymptotic range: 239.0881
## 
## likfit: maximised log-likelihood = -154.5

Valores de AIC e BIC respectivamente para o ml1

aicbic1=c(ml1$AIC,ml1$BIC);aicbic1

## [1] 317.0398 322.5090

Ajuste do modelo “ml2” sem efeito de tendência e com Trasformação.

ml2 <- likfit(INF,lam=0, ini=c(343145.5, 6.66 ),nug= 0)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml2

## likfit: estimated model parameters:
##      beta     tausq   sigmasq       phi 
## " 6.6966" " 0.0000" " 0.1431" "76.3245" 
## Practical Range with cor=0.05 for asymptotic range: 228.6476
## 
## likfit: maximised log-likelihood = -153.8

Valores de AIC e BIC respectivamente para o ml2

aicbic2=c(ml2$AIC,ml2$BIC);aicbic2

## [1] 315.5280 320.9972

Ajuste do modelo “ml3”com efeito de tendência de primeira ordem “1st” e sem trasformação.

ml3 <- likfit(INF, trend="1st", ini=c(343145.5, 6.66 ),nug= 0)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml3

## likfit: estimated model parameters:
##      beta0      beta1      beta2      tausq    sigmasq        phi 
## "1536.547" " -37.656" " -39.559" "   0.000" "1401.524" "   1.705" 
## Practical Range with cor=0.05 for asymptotic range: 5.107397
## 
## likfit: maximised log-likelihood = -139.8

Valores de AIC e BIC respectivamente para o ml3

aicbic3=c(ml3$AIC,ml3$BIC);aicbic3

## [1] 291.5319 299.7357

Ajuste do modelo “ml4”com efeito de tendência de primeira ordem “1st” com trasformação.

ml4 <- likfit(INF, trend="1st",lam=0, ini=c(343145.5, 6.66 ),nug= 0)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml4

## likfit: estimated model parameters:
##     beta0     beta1     beta2     tausq   sigmasq       phi 
## " 7.5220" "-0.0426" "-0.0487" " 0.0000" " 0.0015" " 1.0496" 
## Practical Range with cor=0.05 for asymptotic range: 3.144218
## 
## likfit: maximised log-likelihood = -138

Valores de AIC e BIC respectivamente para o ml4

aicbic4=c(ml4$AIC,ml4$BIC);aicbic4

## [1] 288.0426 296.2464

Ajuste do modelo “ml5” com efeito de tendência na coordenada x sem trasformação.

ml5 <- likfit(INF,trend = ~coords[,1], ini=c(343145.5, 6.66 ),nug= 0)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml5

## likfit: estimated model parameters:
##      beta0      beta1      tausq    sigmasq        phi 
## " 1336.50" "  -44.75" "    0.00" "19806.10" "   22.26" 
## Practical Range with cor=0.05 for asymptotic range: 66.69733
## 
## likfit: maximised log-likelihood = -148.1

Valores de AIC e BIC respectivamente para o ml5

aicbic5=c(ml5$AIC,ml5$BIC);aicbic5

## [1] 306.2378 313.0743

Ajuste do modelo “ml6” com efeito de tendência na coordenada x com trasformação.

ml6 <- likfit(INF,trend = ~coords[,1],lam=0, ini=c(343145.5, 6.66 ),nug= 0)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml6

## likfit: estimated model parameters:
##     beta0     beta1     tausq   sigmasq       phi 
## " 7.2535" "-0.0523" " 0.0000" " 0.0284" "21.4481" 
## Practical Range with cor=0.05 for asymptotic range: 64.25264
## 
## likfit: maximised log-likelihood = -147.8

Valores de AIC e BIC respectivamente para o ml6

aicbic6=c(ml6$AIC,ml6$BIC);aicbic6

## [1] 305.5666 312.4031

Ajuste do modelo “ml7” com efeito de tendência na coordenada y sem trasformação.

ml7 <- likfit(INF,trend = ~coords[,2], ini=c(343145.5, 6.66 ),nug= 0)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml7

## likfit: estimated model parameters:
##      beta0      beta1      tausq    sigmasq        phi 
## " 1194.34" "  -46.80" "    0.00" "38730.95" "   40.83" 
## Practical Range with cor=0.05 for asymptotic range: 122.3032
## 
## likfit: maximised log-likelihood = -149.6

Valores de AIC e BIC respectivamente para o ml7

aicbic7=c(ml7$AIC,ml7$BIC);aicbic7

## [1] 309.1024 315.9388

Ajuste do modelo “ml8” com efeito de tendência na coordenada y com trasformação.

ml8 <- likfit(INF,trend = ~coords[,2],lam=0, ini=c(343145.5, 6.66 ),nug= 0)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml8

## likfit: estimated model parameters:
##     beta0     beta1     tausq   sigmasq       phi 
## " 7.0961" "-0.0565" " 0.0000" " 0.0502" "36.5207" 
## Practical Range with cor=0.05 for asymptotic range: 109.4063
## 
## likfit: maximised log-likelihood = -148.7

Valores de AIC e BIC respectivamente para o ml8

aicbic8=c(ml8$AIC,ml8$BIC);aicbic8

## [1] 307.4990 314.3354

Função matern

Ajuste do modelo inicial “ml9” sem efeito de tendência e sem trasformação.

ml9 <- likfit(INF, ini=c(343145.5 , 2.54 ),nug= 0, kappa = 9.86)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml9

## likfit: estimated model parameters:
##        beta       tausq     sigmasq         phi 
## "  783.652" "  564.608" "50871.521" "    1.305" 
## Practical Range with cor=0.05 for asymptotic range: 14.53213
## 
## likfit: maximised log-likelihood = -150.7

Valores de AIC e BIC respectivamente para o ml9

aicbic9=c(ml9$AIC,ml9$BIC);aicbic9

## [1] 309.4885 314.9577

Ajuste do modelo “ml10” sem efeito de tendência e com Trasformação.

ml10 <- likfit(INF,lam=0, ini=c(343145.5 , 2.54 ),nug= 0, kappa = 9.86)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml10

## likfit: estimated model parameters:
##     beta    tausq  sigmasq      phi 
## "6.4173" "0.0011" "0.3442" "2.4725" 
## Practical Range with cor=0.05 for asymptotic range: 27.52899
## 
## likfit: maximised log-likelihood = -152.2

Valores de AIC e BIC respectivamente para o ml10

aicbic10=c(ml10$AIC,ml10$BIC);aicbic10

## [1] 312.3285 317.7977

Ajuste do modelo “ml11”com efeito de tendência de primeira ordem “1st” e sem trasformação.

ml11 <- likfit(INF, trend="1st", ini=c(343145.5 , 2.54 ),nug= 0, kappa = 9.86)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml11

## likfit: estimated model parameters:
##     beta0     beta1     beta2     tausq   sigmasq       phi 
## "1515.77" " -35.89" " -38.77" "1561.09" "   0.00" "   0.00" 
## Practical Range with cor=0.05 for asymptotic range: 0.0001159668
## 
## likfit: maximised log-likelihood = -147.8

Valores de AIC e BIC respectivamente para o ml11

aicbic11=c(ml11$AIC,ml11$BIC);aicbic11

## [1] 307.5395 315.7432

Ajuste do modelo “ml12”com efeito de tendência de primeira ordem “1st” com trasformação.

ml12 <- likfit(INF, trend="1st",lam=0, ini=c(343145.5 , 2.54 ),nug= 0, 
               kappa = 9.86)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml12

## likfit: estimated model parameters:
##     beta0     beta1     beta2     tausq   sigmasq       phi 
## " 7.5057" "-0.0401" "-0.0489" " 0.0017" " 0.0000" " 0.0000" 
## Practical Range with cor=0.05 for asymptotic range: 0.0001159668
## 
## likfit: maximised log-likelihood = -142.7

Valores de AIC e BIC respectivamente para o ml12

aicbic12=c(ml12$AIC,ml12$BIC);aicbic12

## [1] 297.4945 305.6983

Ajuste do modelo “ml13” com efeito de tendência na coordenada x sem trasformação.

ml13 <- likfit(INF,trend = ~coords[,1], ini=c(343145.5 , 2.54 ),nug= 0,
               kappa = 9.86)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml13

## likfit: estimated model parameters:
##        beta0        beta1        tausq      sigmasq          phi 
## " 1306.8927" "  -44.6548" "  555.9753" "12223.2661" "    0.9811" 
## Practical Range with cor=0.05 for asymptotic range: 10.92419
## 
## likfit: maximised log-likelihood = -147.5

Valores de AIC e BIC respectivamente para o ml13

aicbic13=c(ml13$AIC,ml13$BIC);aicbic13

## [1] 304.9874 311.8239

Ajuste do modelo “ml14” com efeito de tendência na coordenada x com trasformação.

ml14 <- likfit(INF,trend = ~coords[,1],lam=0, ini=c(343145.5 , 2.54 ),nug= 0,
               kappa = 9.86)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml14

## likfit: estimated model parameters:
##     beta0     beta1     tausq   sigmasq       phi 
## " 7.2357" "-0.0533" " 0.0008" " 0.0181" " 1.0026" 
## Practical Range with cor=0.05 for asymptotic range: 11.16285
## 
## likfit: maximised log-likelihood = -147

Valores de AIC e BIC respectivamente para o ml14

aicbic14=c(ml14$AIC,ml14$BIC);aicbic14

## [1] 303.9325 310.7690

Ajuste do modelo “ml15” com efeito de tendência na coordenada y sem trasformação.

ml15 <- likfit(INF,trend = ~coords[,2], ini=c(343145.5 , 2.54 ),nug= 0,
               kappa = 9.86)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml15

## likfit: estimated model parameters:
##        beta0        beta1        tausq      sigmasq          phi 
## " 1140.5829" "  -45.2068" "  523.3132" "19657.4201" "    0.9958" 
## Practical Range with cor=0.05 for asymptotic range: 11.0876
## 
## likfit: maximised log-likelihood = -149.1

Valores de AIC e BIC respectivamente para o ml15

aicbic15=c(ml15$AIC,ml15$BIC);aicbic15

## [1] 308.1299 314.9664

Ajuste do modelo “ml16” com efeito de tendência na coordenada y com trasformação.

ml16 <- likfit(INF,trend = ~coords[,2],lam=0, ini=c(343145.5 , 2.54 ),nug= 0,
               kappa = 9.86)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml16

## likfit: estimated model parameters:
##     beta0     beta1     tausq   sigmasq       phi 
## " 7.0367" "-0.0541" " 0.0008" " 0.0281" " 1.0085" 
## Practical Range with cor=0.05 for asymptotic range: 11.22885
## 
## likfit: maximised log-likelihood = -148.6

Valores de AIC e BIC respectivamente para o ml16

aicbic16=c(ml16$AIC,ml16$BIC);aicbic16

## [1] 307.2130 314.0495

TABELA: Mostra os valores de AIC e BIC dos betas, do tau sigma phi, para os respectivos modelos. Tem-se como base a aceitação do modelo pelo critério de menor AIC ou BIC.

Modelos	AIC	BIC	Beta0	Beta1	Beta2	Tau	Sigma	Phi
ml1	317.0398	322.5090	865.85	—-	—-	0.00	103420.24	79.81
ml2	315.5280	320.9972	6.6966	—-	—-	0.00	0.1431	76.32
ml3	291.5319	299.7357	1536.55	-37.656	-39.559	0.00	1401.524	1.705
ml4	288.0426	296.2464	7.5220	-0.0426	-0.0487	0.00	0.0015	1.0496
ml5	306.2378	313.0743	1336.50	-44.75	—-	0.00	19806.10	22.26
ml6	305.5566	312.4031	7.2535	-0.0523	—-	0.00	0.0284	21.4481
ml7	309.1024	315.9388	1194.34	-46.80	—-	0.00	38730.95	40.83
ml8	307.4990	314.3354	7.0961	-0.0565	—-	0.00	0.0502	36.5207
ml9	309.4885	314.9577	783.652	—-	—-	0.00	50871.521	1.305
ml10	312.3285	317.7977	6.4173	—-	—-	0.00	0.3442	2.4725
ml11	307.5395	315.7432	1515.77	-35.89	-38.77	1561.09	0.00	0.00
ml12	297.4945	305.6983	7.5057	-0.0401	-0.0489	0.0017	0.00	0.00
ml13	304.9874	311.8239	1306.89	-44.655	—-	555.975	12223.266	0.9811
ml14	303.9325	310.7690	7.2357	-0.0533	—-	0.0008	0.0181	1.0026
ml15	308.1299	314.9664	1140.59	-45.207	—-	523.31	19657.32	19657.4
ml16	307.2131	314.0495	7.0367	-0.0541	—-	0.0008	0.0281	1.0085

Melhor modelo

O melhor modelo escolhido foi por meio do critério do menor AIC, portanto o modelo que apresentou esta caracteristica foi o ml3.

Ajuste do modelo

Ajustes dos modelos com as funções de Semivariância.

##############
##   WAVE   ##
##############


salvar=variog(INF)

## variog: computing omnidirectional variogram

plot(salvar,ylab=expression("Semivariância ("*Wave*")"),
     xlab="Distancia (m)", main = "Semivariograma")
lines.variomodel(cov.model="wave", cov.pars=c(343145.5, 6.66 ),nug= 0,kappa =.5)

################
##   MATERN   ##
################

salvar=variog(INF)

## variog: computing omnidirectional variogram

plot(salvar,ylab=expression("Semivariância ("*matern*")"),
     xlab="Distancia (m)", main = "Semivariograma")
lines.variomodel(cov.model="matern", cov.pars=c(343145.5 , 2.54 ),nug= 0,
               kappa = 9.86)

Definindo as bordas arbitrárias

Bordas do modelo.

points(INF, pt.div="quartiles", xlab="Leste", ylab="Norte")

setwd("F:/Estatistica Espacial")
borda=read.table("bordaINF.txt",h=T)
INF$borders <- with(borda, cbind(x,y))

points(INF, pt.div="quartiles", xlab="leste", ylab="norte")

plot(INF,lowess = T)

Tomando apenas como base o gráfico acima, pode-se observar que os dados não apresentam normalidade.

Boxcox

Por meio do Boxcox verifica-se se o 1 está contido no intervalo, se caso estiver,conclui-se que os dados nececitam de trasformação.

boxcox(INF,lam=seq(-3,2,1/10))

Interpolação

par(mfrow=c(1,1))
ml3 <- likfit(INF, trend="1st", ini=c(343145.5, 6.66 ),nug= 0)

## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.

ml3

## likfit: estimated model parameters:
##      beta0      beta1      beta2      tausq    sigmasq        phi 
## "1536.547" " -37.656" " -39.559" "   0.000" "1401.524" "   1.705" 
## Practical Range with cor=0.05 for asymptotic range: 5.107397
## 
## likfit: maximised log-likelihood = -139.8

gr=pred_grid(INF$borders,by=0.2)
points(INF)
points(gr,col=2,pch=19,cex=.3)

gr0=gr[.geoR_inout(gr,INF$borders),]
points(gr0,col=3,pch=19,cex=0.3)

Processo de Krigagem

O processo de Krigagem permite construir mapas de incerteza espacial, no qual fornece uma visão global da congruência dos resultados.

kc <- krige.conv(INF, loc=gr, bor=borda, krige=krige.control(obj=ml3))

## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with mean given by a 1st order polynomial on the coordinates
## krige.conv: Kriging performed using global neighbourhood

names(kc)

## [1] "predict"      "krige.var"    "beta.est"     "distribution"
## [5] "message"      "call"

kc$pred[1:10]

##  [1] 1043.831 1036.046 1028.895 1025.360 1023.141 1020.784 1206.516
##  [8] 1197.951 1050.759 1043.152

image(kc,main="Mapa de Interpolação")

image(kc, main="Mapa de Interpolação",col=gray(seq(1,0.2,len=10)), x.leg=c(-220,-50), y.leg=c(150,180))

kc1 <- krige.conv(INF, loc=gr, bor=borda, krige=krige.control(obj=ml3))

## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with mean given by a 1st order polynomial on the coordinates
## krige.conv: Kriging performed using global neighbourhood

plot(kc$pred, INF$pred, asp=10)
abline(0,1, col=2)

image(kc1, main="Mapa de Interpolação",col=gray(seq(1,0.2,len=20)), x.leg=c(-220,-50), y.leg=c(150,180))

contour(kc1,main="Mapa de Interpolação",f=T,col= terrain.colors(20), nlev = 20)

numlevel=10
contour(kc1,main="Mapa de Interpolação",c=T,col= terrain.colors(50), nlev = 50)

numlevel=10