Faire de la géostatistique avec R

Comme vous l’avez vu en cours nous allons utiliser le package geoR pour faire notre analyse spatiale. Nous allons reprendre l’analyse du jeu de donnees de concentration en Calcium vu en cours. Attention toutefois je ne commenterai pas les résultats car nous l’avons déjà fait en cours. Le but de ce document est simplement de vous aider dans l’utilisation de geoR (et autres)…

1. Analyse descriptive

Après avoir importé les données et vérifié qu’il n’y avait pas de valeurs abérantes, nous faisons donc un symbol plot

library(geoR)
data(ca20)
## Mon code pour le symbol plot
par(ps = 14, mar = c(4, 5, 0.5, 0.5))
plot(ca20$borders, type = "l", xlab = "Longitude", ylab = "Latitude",
     asp = 1)
lines(ca20$reg1, lty = 2)

lines(ca20$reg2, lty = 2)
lines(ca20$reg3, lty = 2)

symbols(ca20$coord, circles = abs(ca20$data-mean(ca20$data)),
        bg = c("blue", "red")[1 + as.numeric(ca20$data > mean(ca20$data))],
        inches = FALSE, add = TRUE)
text(ca20$coord, labels = round(ca20$data, 1), cex = 0.5)

par(mfrow = c(1, 4), ps = 20, mar = c(4, 5, 0.75, 0.5))

plot(ca20$coord[,1], ca20$data, xlab = "Longitude", ylab = "Ca")
lines(lowess(ca20$coord[,1], ca20$data), col = "orange")

plot(ca20$coord[,2], ca20$data, xlab = "Latitude", ylab = "Ca")
lines(lowess(ca20$coord[,2], ca20$data), col = "orange")

plot(ca20$covariate[,"altitude"], ca20$data, xlab = "Altitude", ylab = "Ca")
lines(lowess(ca20$covariate[,"altitude"], ca20$data), col = "orange")

plot(ca20$covariate[,"area"], ca20$data, xlab = "Sous-region", ylab = "Ca")

Cela dit puisque le jeu de données Ca20 est au format geodata, il existe une méthode pour cela (mais en pratique ce ne sera pas le cas à moins de convertir vous même vos données au format geodata)

plot(ca20)

Passons au variogramme empirique.

vario.emp <- variog(coords = ca20$coords, data = ca20$data)

variog: computing omnidirectional variogram

plot(vario.emp)

## On aurait pu faire directement plot(variog(ca20)) car ca20 est de classe geodata

Pour avoir un variogramme empirique directionnel

plot(variog4(coords = ca20$coords, data = ca20$data))

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

2. Modélisation paramétrique

Deux approches sont alors possibles : soit l’on ajuste un variograme paramétrique via les moindres carrés (généralisés ou pas) soit en ajustant directement un processus gaussien.

fit.vario <- variofit(vario.emp, cov.model = "stable", fix.nugget = TRUE)##ici on fait confiance au package pour qu'il trouve de bonnes valeurs initiales, sinon on fait comme dit en cours ;-)

variofit: covariance model used is powered.exponential 
variofit: weights used: npairs 
variofit: minimisation function used: optim 

initial values not provided - running the default search

variofit: searching for best initial value ... selected values:
              sigmasq  phi     tausq kappa
initial.value "153.65" "175.1" "0"   "0.5"
status        "est"    "est"   "fix" "fix"
loss value: 3360490.37531911 

plot(vario.emp, ylim = c(0, 200))
lines(fit.vario, col = "seagreen3")

Bon d’accord mais c’est un peu nul ce que l’on a fait non car il y a un effet de la sous-région ! On en tient compte alors

vario.emp2 <- variog(ca20, trend = ~ area)

variog: computing omnidirectional variogram

plot(vario.emp2)
fit.vario2 <- variofit(vario.emp2, fix.nugget = TRUE)

variofit: covariance model used is matern 
variofit: weights used: npairs 
variofit: minimisation function used: optim 

initial values not provided - running the default search

variofit: searching for best initial value ... selected values:
              sigmasq  phi     tausq kappa
initial.value "106.36" "175.1" "0"   "0.5"
status        "est"    "est"   "fix" "fix"
loss value: 1762223.86576382 

lines(fit.vario2, col = "seagreen3")

Maintenant si je veux passer par l’ajustement d’un processus gaussien

fit.gp <- likfit(ca20, trend = ~ area, fix.nugget = TRUE, ini.cov.pars = c(100, 100), cov.model = "stable", fix.kappa = FALSE)

---------------------------------------------------------------
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.

Comparons à la fonction que nous avons fait en cours et disponible ici

url <- "http://imag.umontpellier.fr/~ribatet/docs/Geostatistique/fitgp.R"
source(url)
fit.gp.mine <- myfitgp(data = ca20$data, coords = ca20$coords,
                       trend = ~ area, covariables = ca20$covariate, iterlim = 1000)

Je vous laisse voir que nous trouvons les mêmes résultats. Conclusion on est trop fort non ?

3. Krigeage

Le plus simple consiste à passer un variogramme ajusté ou un processus gaussien ajusté directement à la fonction faisant le krigeage.

gr <- pred_grid(ca20$borders, by = 2.5)## points ou l'on veut predire
gr0 <- polygrid(gr, borders = ca20$borders, bound = TRUE)## on se restreint a la region d'etude et on n'ira pas en dehors
ind.reg <- numeric(nrow(gr0))
ind.reg[.geoR_inout(gr0, ca20$reg1)] <- 1
ind.reg[.geoR_inout(gr0, ca20$reg2)] <- 2
ind.reg[.geoR_inout(gr0, ca20$reg3)] <- 3
ind.reg <- as.factor(ind.reg)##ici j'ai affecte le numero de la sous region a chaque point de prediction
KC <- krige.control(obj.model = fit.gp, trend.d = ~ area, trend.l = ~ ind.reg)
ca20pred <- krige.conv(ca20, loc = gr, krige = KC)

krige.conv: results will be returned only for prediction locations inside the borders
krige.conv: model with mean defined by covariates provided by the user
krige.conv: Kriging performed using global neighbourhood 

## Ouf on peut faire le graphique
par(mfrow = c(1, 2), mar = c(4, 5, 0.5, 0.5), ps = 14)
##plot(variog(ca20), xlab = expression(h), ylab = expression(gamma(h)))
##f <- function(h) (16.77 + 135.17) - (16.77 * as.numeric(h == 0) + 135.17) * exp(-h / 159.49)
##curve(f, from = 0, to = 1500, col = "orange", add = TRUE)
image(ca20pred, loc = gr, y.leg = c(4790, 4840), x.leg = c(4930, 5350),
      col = heat.colors(64), breaks = seq(15, 80, length = 65),
      zlim = c(15, 80))
lines(ca20$reg1, lty = 2)

lines(ca20$reg2, lty = 2)
lines(ca20$reg3, lty = 2)

text(ca20$coord, labels = ca20$data, cex = 0.75)
image(ca20pred, loc = gr, y.leg = c(4790, 4840), x.leg = c(4930, 5350),
      val = sqrt(ca20pred$krige.var), col = heat.colors(64), breaks = seq(0, 13, length = 65),
      zlim = c(0, 13))

lines(ca20$reg1, lty = 2)
lines(ca20$reg2, lty = 2)

lines(ca20$reg3, lty = 2)
points(ca20$coord, pch = 15)