Pacotes utilizados e leitura do banco de dados

library(sp)
library(geoR) 
library(gstat)
library(tidyverse)
library(sf)
library(maptools)
library(spdep)
library(RColorBrewer)
library(classInt)
library(gridExtra)
library(nortest)

dados_sim.prova2 <- read.csv("/cloud/project/Prova2Spatial/dados_sim-prova2.txt", sep="")
#View(dados_sim.prova2)

Banco de dados

O banco de dados referente a seunda avaliação da disciplina, apresenta 93 observações sobre a resistência do solo à penetração na camada de \(0 - 0.10 cm\), medidas em MPa. É suposto que os dados foram coletados em uma área de 167 hectares, onde o clima apresenta como mesotérmico e super úmido, temperaturas moderadas, chuvas bem distribuídas e verão quente. O solo é classificado como argiloso.

Análise descritiva Espacial

Nesta seção será feita uma análise descritiva espacial, no qual o principal objetivo é de vizualizar e entender a natureza dos dados apresentados. Esse estudo pode ser feito tanto olhando apena para a variável, como também espacialmente falando, já que no banco de dados temos essa variável georeferenciada, ou seja, com as coordenada “x” e “y”.

Análise descritiva usual

Antes de se utilizar das vizualizações dos mapas, é preciso utilizar uma análise descritiva da variável, afim de olhar apenas para as medições e não suas posições geográficamente. No boxplot, podemos afirmar que existe a presença de 2 outliers, ou seja, medições que tiveram uma resistência alta. No histograma, é mais notável a questão da assimetria, no qual os dados apresentam uma leve assimetria positiva, ou seja, uma maior concentração para valores menores dessas medições de resistência.

data<-dados_sim.prova2$data %>% as.data.frame()
colnames(data)<- paste(c("Res"))

p1<-data %>% ggplot() +
geom_boxplot(aes(x="", y = Res), fill = "grey", color = "darkblue") +
labs(y="Restistência")


p2<-data %>%  ggplot() +
geom_histogram(aes(x = Res), fill = "lightblue", color = "darkblue", binwidth = 0.7) +
labs(x="Restistência", y="Frequência")

grid.arrange(p1,p2, ncol=2,nrow=1)

Trabalhando com dados georeferenciados

Como estamos agora utilizando as coordenadas a nosso favor para estudar o comportamento da variável no mapa, temos que mudar a forma que utilizamos o banco de dados. Originalmente foi utilizado ele como um data frame, porém, quando estamos falando em geoestatísica, é preciso tratar esses dados de forma diferenciada. Essa forma diferenciada seria em transforma esse data frame em uma outra classe, denominada de “SpatialPointsDataFrame”.

Após tranformar o data.frame original, como os dados apresentados são da turquia, foi feita uma busca na internet e foi encontrado um sistem de referência “epsg:5636”, nele foi usado para plotar a figura e vizualizar a distribuição da variável em relação as suas coordenadas.

#y -> Latitude
#x -> Longitude
coordinates(dados_sim.prova2)<-~x+y
proj4string(dados_sim.prova2) <-  CRS("+init=epsg:5636")
plot(dados_sim.prova2, axes=TRUE, pch = 20, col = "darkblue", asp = 1)
grid(lty = 1, col = "darkgray")

#dados_sim.prova2@coords

Vale notar no gráfico acima que os dados estão no sistema de coordenada em UTM, outra possibilidade é de usarmos o sistema em formato “longlat”. Utilizando os códigos abaixo e vizualizando, os dados estão presentes em um terreno com 39.6838° à 40.9938° Leste e 37.1° à 37.6601° Norte.

#coordinates(dados_sim.prova2)<-~y+x
proj4string(dados_sim.prova2) <-  CRS("+proj=longlat")
plot(dados_sim.prova2, axes=TRUE, pch = 20, col = "darkblue", asp = 1)
grid(lty = 1, col = "darkgray")

bbox(dados_sim.prova2@coords)
##       min     max
## x 39.6838 40.9938
## y 37.0731 37.6601

Porém, nesses dois gráficos acima, foi visto apenas as coordenadas da variável, não incluindo qualquer quantidade ou variabilidade da mesma.

Nesse gráfico, foi introduzido a variável em estudo, nela podemos ver que existe uma maior resistência ao solo na parte nordeste e noroeste do mapa, enquanto no centro da região, essa resistência no solo à penetração diminui. Essa característica pode se referir pelas característica geográficas da região, onde podemos supor que na região nordesete e noroste refere-se a lugares mais montanhosos, enquanto no centro é mais plano.

plotvar <- dados_sim.prova2$data
nclr <- 4
plotclr <- brewer.pal(nclr,"YlOrRd")
class <- classIntervals(plotvar, nclr, style="quantile")
colcode <- findColours(class, plotclr)
plot(dados_sim.prova2, axes=TRUE)
points(dados_sim.prova2@coords[,1], dados_sim.prova2@coords[,2], pch=16, col=colcode, cex=2)
points(dados_sim.prova2@coords[,1], dados_sim.prova2@coords[,2], cex=2)
legend("topright", legend=names(attr(colcode, "table")),
    fill=attr(colcode, "palette"), cex=0.6, bty="n")

Estudo da anisotropia.

É uma característica que representa mudança daquele fênomeno em estudo, em relação as direção em que são medidas. A anisotropia implica que a dependência espacial não é uma função apenas das distâncias entre as observações, mas sim em que direção elas são observadas espacialmentem nos locais de amostragem. A anisotropia pode ser classficada em três tipos:

Com isso, precisamos construir o gráfico de semivariâncias, mais conhecido como o semivariograma. Foi utilizado abaixo o pacote geoR. O cutoff utilizado foi de \(0.8\). Podemos perceber nesse semivariograma, que utilizando a distãncia máxima de 0.8, os dados apresentam uma características bem linear até certo ponto.

dados_sim<-as.data.frame(dados_sim.prova2)
dados_sim2<-as.geodata(obj = dados_sim, coords.col = 1:2, data.col = 3)
v1<- variog(dados_sim2, coords = dados_sim2$coords, data = dados_sim2$data, estimator.type = "classical", max.dist=0.8, pairs.min = 30)
## variog: computing omnidirectional variogram
plot(v1, main="Variograma Empírico dos dados", ylab="Semivariância", xlab="Distâncias (Metros)")

Agora, em respeito a anisotropia, precisamos saber se existe mudança significativa nesse semivariograma quando olhamos em diferentes direções. Utilizando a função variog4 do pacote geoR, é obtido o gráfico abaixo.

v2<- variog4(dados_sim2, coords = dados_sim2$coords, data = dados_sim2$data, estimator.type = "classical", max.dist=0.8, pairs.min = 30, direction = c(0, pi/4, pi/2, 3*pi/4))
## variog: computing variogram for direction = 0 degrees (0 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogram
plot(v2, main="Variograma Empírico dos dados para diferentes direções", ylab="Semivariância", xlab="Distâncias (Metros)")

Neste caso, podemos observar que nas direções 0°, 45° e 90° apresentam um alcance menor que olhando para a direção de 90°, sendo a de 0° a de menor alcance. Em relação a contribuição, as direções que apresentam menor valor foram a de 0° e 90°. A direção de 45° foi a que houve maior valor na contribuição.

De forma geral, podemos dizer que os dados apresentam uma anisotropia mista, apresentando tanto mudanças na contribuicões e nos alcances, em relação as direções.

Distribuição dos dados.

Nesta seção será feito o estudo, para encontrar a distribuição de probabilidade que mais caracteriza o comportamento dos dados, em que, isso é de extrema importância para problemas posteriores.

A primeira suposição a se fazer é da normalidade dos dados, assim, foi feito

#Avaliando normalidade dos dados
# Gráficamente
attach(data)
par(mfrow=c(2,2))
hist(Res, main="Histograma", col="grey", data=data, xlab = "Resistência", ylab="Frequência")
boxplot(Res, main="Boxplot", col="grey", data=data, ylab = "Resistência")
qqnorm(Res, main="Gráfico de Quantís", ylab="Quantís Amostrais ", xlab="Quantís Teóricos")
qqline(Res, lty = 2, col = "red")
plot(density(Res), main="Gráfico da Densidade",  data=data, ylab="Densidade")

Como a análise apenas por gráficos pode ser um pouco subjetiva ou tendenciosa, então a maneira mais prática é de utilizarmos teste estatísticos de aderência. Nesse caso, foi feito cerca de seis testes, entre eles o mais conhecido pela literatura, o teste de Kolmogorov-Smirnov. O critério de decisão é de escolher o maior fator em comum desses testes, ou seja, escolher a decisão apontada pela maioria.

#Estudo de normalidade da variável Resposta
t1 <- ks.test(Res, "pnorm", mean(Res), sd(Res))  # KS 
t2 <- lillie.test(Res)  # Lilliefors
t3 <- cvm.test(Res)  # Cramér-von Mises
t4 <- shapiro.test(Res) # Shapiro-Wilk
t5 <- sf.test(Res) # Shapiro-Francia
t6 <- ad.test(Res) # Anderson-Darling

# Tabela de resultados
testes <- c(t1$method, t2$method, t3$method, t4$method, t5$method, t6$method) 
estt <- as.numeric(c(t1$statistic, t2$statistic, t3$statistic, t4$statistic, t5$statistic, t6$statistic))
valorp <- c(t1$p.value, t2$p.value, t3$p.value, t4$p.value, t5$p.value, t6$p.value)
resultados <- cbind(estt, valorp)
rownames(resultados) <- testes
colnames(resultados) <- c("Estatística", "P-valor")
print(round(resultados,5), digits = 5)
##                                                Estatística P-valor
## One-sample Kolmogorov-Smirnov test                 0.10536 0.23603
## Lilliefors (Kolmogorov-Smirnov) normality test     0.10536 0.01265
## Cramer-von Mises normality test                    0.27187 0.00068
## Shapiro-Wilk normality test                        0.92957 0.00009
## Shapiro-Francia normality test                     0.93221 0.00027
## Anderson-Darling normality test                    1.86280 0.00009

Com isso, temos a evidência de rejeitar a hipótese nula que os dados possuem distribuição normal, para um nível de significância \(\alpha=5\%\)

Porém podemos utilizar de uma transformação de variável, uma tranformação válida é a logaritma, com ela o comportamento dos dados melhora em relação a questão da normalidade. Observando os gráficos abaixo, pelo histograma, os dados apararentemente estão se comportando de forma simétrica, com algumas excessões para alguns valores menores do MPa, levando a ideia de estarmos nos deparando com distribuições de caudas mais pesadas. Ainda com a possibilidade da normalidade, temo os boxplot que comporta de maneira simétrica, com a mediana não viesada para nenhum dos quartis. Um dos gráficos mais interessantes a se falar quando estamos querendo estudar a respeito da normalidade dos dados, é o gráfico de quantis, no qual, ele compara os quantis da variável em estudo com os quantis de uma distribuição normaL, em que, se os dados permanecerem alinhados a linha vermelha tracejada , maior a evidência (de forma vizual) que os dados são normalmente distribuídos. Neste mesmo gráfico, podemos perceber que houve um bom comportamento em relação a linha vermelha, com excessões aos valores mais abaixo.

#Avaliando normalidade dos dados
# Gráficamente
data$Res_new<- log(Res)
par(mfrow=c(2,2))
attach(data)
hist(Res_new, main="Histograma", col="grey", data=data, xlab = "Resistência", ylab="Frequência")
boxplot(Res_new, main="Boxplot", col="grey", data=data, ylab = "Resistência")
qqnorm(Res_new, main="Gráfico de Quantís", ylab="Quantís Amostrais ", xlab="Quantís Teóricos")
qqline(Res_new, lty = 2, col = "red")
plot(density(Res_new), main="Gráfico da Densidade",  data=data, ylab="Densidade")

Pra confirmar a análise acima

#Estudo de normalidade da variável Resposta
t1 <- ks.test(Res_new, "pnorm", mean(Res_new), sd(Res_new))  # KS 
t2 <- lillie.test(Res_new)  # Lilliefors
t3 <- cvm.test(Res_new)  # Cramér-von Mises
t4 <- shapiro.test(Res_new) # Shapiro-Wilk
t5 <- sf.test(Res_new) # Shapiro-Francia
t6 <- ad.test(Res_new) # Anderson-Darling

# Tabela de resultados
testes <- c(t1$method, t2$method, t3$method, t4$method, t5$method, t6$method) 
estt <- as.numeric(c(t1$statistic, t2$statistic, t3$statistic, t4$statistic, t5$statistic, t6$statistic))
valorp <- c(t1$p.value, t2$p.value, t3$p.value, t4$p.value, t5$p.value, t6$p.value)
resultados <- cbind(estt, valorp)
rownames(resultados) <- testes
colnames(resultados) <- c("Estatística", "P-valor")
print(round(resultados,5), digits = 5)
##                                                Estatística P-valor
## One-sample Kolmogorov-Smirnov test                 0.05769 0.89858
## Lilliefors (Kolmogorov-Smirnov) normality test     0.05769 0.62585
## Cramer-von Mises normality test                    0.05742 0.40529
## Shapiro-Wilk normality test                        0.97890 0.13722
## Shapiro-Francia normality test                     0.98293 0.22707
## Anderson-Darling normality test                    0.47477 0.23528

Utilizando agora os dados transformados para a obtenção do variograma.

dados_sim<-as.data.frame(dados_sim.prova2)
dados_sim$data<-log(dados_sim$data)
dados_sim2<-as.geodata(obj = dados_sim, coords.col = 1:2, data.col = 3)
v1_new<- variog(dados_sim2, coords = dados_sim2$coords, data = dados_sim2$data, estimator.type = "classical", max.dist=0.8, pairs.min = 30)
## variog: computing omnidirectional variogram
plot(v1_new, main="Variograma Empírico dos dados transformados", ylab="Semivariância", xlab="Distâncias (Metros)")

Escolha do melhor modelo

Nesta seção será estudado qual modelo mais se ajusta aos dados do semivariograma acima. Vai ser testa cerca de 4 modelo para cada tipo de estimador, a análise vai ser feita vizualmente e por outros critérios como o efeito pepita relativo e técnicas de validação cruzada.

Na função variofit os modelos de variogramas podem ser ajustados usando mínimos quadrados ponderados ou ordinários.

m1_wls<-variofit(v1_new, cov.model = "spherical")
m2_wls<-variofit(v1_new, cov.model = "gaussian")
m3_wls<-variofit(v1_new, cov.model = "exponential")
m4_wls<-variofit(v1_new, cov.model = "matern", kappa = 0.30)

m1_ols<-variofit(v1_new, cov.model = "spherical", weights = "equal")
m2_ols<-variofit(v1_new, cov.model = "gaussian", weights = "equal")
m3_ols<-variofit(v1_new, cov.model = "exponential", weights = "equal")
m4_ols<-variofit(v1_new, cov.model = "matern", kappa = 0.30 , weights = "equal")
## variofit: covariance model used is spherical 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq  kappa
## initial.value "0.1"   "0.61" "0.01" "0.5"
## status        "est"   "est"  "est"  "fix"
## loss value: 0.343426539785544
## variofit: covariance model used is gaussian 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq  kappa
## initial.value "0.1"   "0.37" "0.03" "0.5"
## status        "est"   "est"  "est"  "fix"
## loss value: 0.413981883301151
## variofit: covariance model used is exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "0.13"  "0.37" "0"   "0.5"
## status        "est"   "est"  "est" "fix"
## loss value: 0.503958926006976
## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq  kappa
## initial.value "0.1"   "0.12" "0.01" "1.5"
## status        "est"   "est"  "est"  "fix"
## loss value: 0.417196260032547
## variofit: covariance model used is spherical 
## variofit: weights used: equal 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq  kappa
## initial.value "0.1"   "0.61" "0.01" "0.5"
## status        "est"   "est"  "est"  "fix"
## loss value: 0.00130100222039168
## variofit: covariance model used is gaussian 
## variofit: weights used: equal 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq  kappa
## initial.value "0.1"   "0.25" "0.01" "0.5"
## status        "est"   "est"  "est"  "fix"
## loss value: 0.00145721289611095
## variofit: covariance model used is exponential 
## variofit: weights used: equal 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq kappa
## initial.value "0.13"  "0.37" "0"   "0.5"
## status        "est"   "est"  "est" "fix"
## loss value: 0.00171059642369551
## variofit: covariance model used is matern 
## variofit: weights used: equal 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi    tausq  kappa
## initial.value "0.1"   "0.12" "0.01" "1.5"
## status        "est"   "est"  "est"  "fix"
## loss value: 0.0015611615032201

Fazendo agora os plot, cada plot vai ter o semivariograma e 4 modelos correspondentes.

Para os modelos utilizando o estimador de mínimos quadrados ponderados, vizualmente os que tiveram melhor resultado foram o “Esférico” e o “Gaussiano”

Já para os modelos utilizando o estimador de mínimos quadrados ordinários, vizualmente os que tiveram melhor resultado foram o “Esférico” e o “Gaussiano”.

Com isso, selecionando esses dois modelos e fazendo validação cruzada pra esses dois modelos.

Os resultados gráficos são mostrados para os resultados da validação cruzada em que foi usada a estratégia de deixar de fora, combinada com as estimativas wls e ols para os parâmetros. Os resíduos da validação cruzada são obtidos subtraindo os dados observados menos o valor previsto. Os resíduos padronizados são obtidos dividindo-se pela raiz quadrada da variação de previsão (‘variação de krigagem’).

Comparando apenas os modelos com os estimadores de mínimos quadrados ponderados, temos que o modelo esférico se adequou mais aos dados, tanto pelo histograma das predições como pelos resíduos. O modelo gaussiano teve bons resultados, porém, apresentou mais valores concentrados nos gráficos subjacentes.

#m1_wls -> Esférico
#m2_wls -> Gaussiano
m1.xv.wls <- xvalid(dados_sim2, model = m1_wls)
## xvalid: number of data locations       = 93
## xvalid: number of validation locations = 93
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 
## xvalid: end of cross-validation
m2.xv.wls <- xvalid(dados_sim2, model = m2_wls)
## xvalid: number of data locations       = 93
## xvalid: number of validation locations = 93
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 
## xvalid: end of cross-validation
##
## Plotting results
##
par.ori <- par(no.readonly = TRUE)
##
par(mfcol=c(5,2), mar=c(2.3,2.3,.5,.5), mgp=c(1.3, .6, 0))
plot(m1.xv.wls, main = "Modelos WLS Esférico")

par(mfcol=c(5,2))
plot(m2.xv.wls, main = "Modelos WLS Gaussiano")

##
par(par.ori)
#

Observando agora os modelos utilizando os métodos de mínimos quadrados ordinários, temos que o modelo esférico apresenta um histograma mais simétrico e analogamente à comparação anterior, os outros gráficos ficaram mais comportados pelo modelo esférico.

#m1_ols -> Esférico
#m2_ols -> Gaussiano
m1.xv.ols <- xvalid(dados_sim2, model = m1_ols)
## xvalid: number of data locations       = 93
## xvalid: number of validation locations = 93
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 
## xvalid: end of cross-validation
m2.xv.ols <- xvalid(dados_sim2, model = m2_ols)
## xvalid: number of data locations       = 93
## xvalid: number of validation locations = 93
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 
## xvalid: end of cross-validation
##
## Plotting results
##
par.ori <- par(no.readonly = TRUE)
##
par(mfcol=c(5,2), mar=c(2.3,2.3,.5,.5), mgp=c(1.3, .6, 0))
plot(m1.xv.ols, main = "Modelos OLS Esférico")

par(mfcol=c(5,2))
plot(m2.xv.ols, main = "Modelos OLS Gaussiano")

##
par(par.ori)
#

Interprete os parâmetros estimados considerando o modelo selecionado no item anterior, e comente sobre a variabilidade espacial.

O dois modelos selecionados foram o esférico utilizando o método ols e o esférico usando o wls. Nesses dois modelos os resultados são quase idênticos.

Observando os resultados dos parâmetros dos dois modelos restantes, temos que no ols um parâmetro \(\phi_2 = 0.6583\) sendo a contribuição. O efeito \(\phi_1 = 0\), denominado de efeito pepita. O parâmetro sigmasq que é o alcance, teve um valor de 0.1225. Pela variabilidade dos dados em questão, foi preferível escolher o modelo esférico usando ols por apresentar uma contribuição maior em relação ao outro.

print(m1_ols)
## variofit: model parameters estimated by OLS (ordinary least squares):
## covariance model is: spherical
## parameter estimates:
##   tausq sigmasq     phi 
##  0.0000  0.1225  0.6583 
## Practical Range with cor=0.05 for asymptotic range: 0.6583002
## 
## variofit: minimised sum of squares = 5e-04
print(m1_wls)
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## parameter estimates:
##   tausq sigmasq     phi 
##  0.0000  0.1224  0.6559 
## Practical Range with cor=0.05 for asymptotic range: 0.6559095
## 
## variofit: minimised weighted sum of squares = 0.1378

Com esse modelo, podemos falar agora da sua dependência espacial, uma media interessantes para esse caso é o efeito pepita relativo, dado pela equação abaixo: \[ \mathrm{EPR} = \left(\dfrac{\phi_1}{\phi_1 + \phi_2}\right) \times 100 \] É uma medida que indica o nível grau de dependência do modelo em estudo. A interpretação para esse índice é bem fácil e prática, dado por:

\[ \mathrm{Decisão} = \left \{ \begin{matrix} EPR < 0,5, & \mbox{Forte dependência espacial} \\ 0,5 ≤ EPR <0,75 , & \mbox{ Média dependência espacial} \\ EPR >0,75, & \mbox{Fraca dependência espacial} \\ \end{matrix} \right. \] Como já esperado pelos análises feita anteriormente, o EPR para esse modelo é de 0, devido ao fato de ter \(\phi_1\) (efeito pepita) = 0.

#Cáculo do EPR
phi1<-m1_ols$nugget
phi2<-m1_ols$cov.pars[2]
EPR <-((phi1)/(phi1+ phi2) ) *100 
print(EPR)
## [1] 0

Obtenha os valores preditos por krigagem e apresente o mapa (apresente também o mapa do desvio padrão dos valores preditos). Interprete.

Utilizando o código abaixo, podemos vizualizar as duas imagens referentes a predições em coordenadas não amostradas, utilizando a krigagem. O tipo de krigagem foi a ordinária, type.krige = "ok

Podemos dizer que a respeito das predições, que para lugares mais afastados indicam que que os valores da resistência ao solo em MPa, vai ser mais elevado do que no interior do terreno. A respeito da variabilidade, para as localizações do interior, a variabilidade da predição é bem baixa, para valores extremos, a variabilidade aumenta consideravelmte. Uma dessas causas seria por estar possivelmente em um terreno mais elevado ou algum tipo de solo específico na região, que pode afetar nas medições do solo.

#args(krige.control)
#args(output.control)
pred.grid <- expand.grid(seq(39.4, 41.2, l = 40), seq(37,37.7 , l = 40))

KC <- krige.control(obj.m = m1_ols, trend.d="1st", trend.l="1st", type.krige = "ok")
OC <- output.control(n.pred=1000, quantile=c(0.10, 0.25, 0.5, 0.75, 0.90), threshold = 350, simulations.predictive = 2) 
dados_sim2.kc <- krige.conv(dados_sim2, loc=pred.grid, krige = KC, output = OC)
## krige.conv: model with constant mean
## krige.conv: sampling from the predictive distribution (conditional simulations)
## krige.conv: Kriging performed using global neighbourhood
names(dados_sim2.kc)
##  [1] "predict"                   "krige.var"                
##  [3] "beta.est"                  "distribution"             
##  [5] "simulations"               "mean.simulations"         
##  [7] "variance.simulations"      "quantiles.simulations"    
##  [9] "probabilities.simulations" "sim.means"                
## [11] ".Random.seed"              "message"                  
## [13] "call"
par(mfrow=c(1,2))
image(dados_sim2.kc ,val= dados_sim2.kc$predict ,loc = pred.grid, xlab = "Coord X", ylab = "Coord Y",main="Gráfico para os valores preditos")
image(dados_sim2.kc ,val= sqrt(dados_sim2.kc$krige.var) ,loc = pred.grid, xlab = "Coord X", ylab = "Coord Y", main="Gráfico para a variância da previsões")

Além do que foi mencionado nos itens anteriores, o que mais pode ser feito para analisar os dados apresentados? Apresente os resultados.

Outra opção a ser feita é no caso da krigagem, explorar mais os recursos que a função do pacote geoR nos proporciona. Um desses casos é de utilizar simulação de dados geoestatísticos gaussiana, que no caso dessa variável transformada, apresenta-se distribuição normal.

A simulação 1 apresentou estar bem mais próxima em relação a krigagem ordinária utilizada anteriormente.

#args(krige.control)
#args(output.control)
par(mfrow=c(1,2))
image(dados_sim2.kc ,val= dados_sim2.kc$simulation[,1] ,loc = pred.grid, xlab = "Coord X", ylab = "Coord Y",main="Simulação 1")
image(dados_sim2.kc ,val= dados_sim2.kc$simulation[,2] ,loc = pred.grid, xlab = "Coord X", ylab = "Coord Y", main="Simulação 2")