---
title: "Effet Zéro par Gibbs"
author: "Xavier FREULON, Didier Renard"
date: "20 novembre 2021"
output:
  pdf_document: default
  html_document: default
---

# Introduction

Ce test a pour but de démontrer le procédé de simulation avec des données censurées (effet zéro). Il démontre l'impact négatif lorsque ces données nulles ne sont pas traitées correctement.

```{r setup, include=FALSE}
library(RGeostats)
rm(list = ls())
set.seed(32241)
constant.define("asp",1)
```

# Jeu de données testé

Les données sont générées sur une grille régulière (100 x 100) couvrant un carré de taille 1.

```{r Donnees}
nx <- rep(100,2)
dx <- 1/nx
x0 <- rep(0.0,2)
grd <- db.create(flag.grid = TRUE, nx = nx, dx = dx, x0 = x0)
```


Les valeurs sont simulées suivant un variogramme exponentiel de portée 0.25 et de palier unitaire.

```{r Modele}
m <- model.create(vartype = 2, range = 0.25, sill = 1.0)
grd <- simtub(,dbout = grd, model = m, seed = 1234, 
              nbsimu = 1, nbtuba = 1000)
```

Une sélection de *np* valeurs est effectuée. Les échantillons sélectionnés constituent la nouvelle base de données *data*.

```{r Selection}
np = 1000
idx <- sample(x = 1:grd$nech, size = np, replace = FALSE)
data = db.reduce(grd,ranks=idx)
data <- db.rename(data, names = "Simu.V1.S1", newnames = "Y")
data <- db.locate(data, names = "Y", loctype = "z", flag.locnew = TRUE )
plot(grd,title="Simulation exhaustive sur grille et sélection des données")
plot(data, add = TRUE)
```

# Contraintes

Un nouveau champ de valeurs *Yc* est créé à partir des données (après transformation en gaussiennes par Normal Score). La proportion de données censurées est fixé par la proportion (*prop*), ce qui correspond à la coupure gaussienne *yc*.

```{r Coding_data, echo = TRUE, eval = TRUE}
prop <- 0.50
data <- anam.z2y(db = data, anam = NA, names = "Y", radix = "Y")
Y  <- db.extract(data, "Y.Y", flag.compress = FALSE)
yc <- qnorm(p = prop, mean = 0, sd = 1)
cat("Proportion =",prop,"- Seuil gaussian yc =",yc,"\n")
Yc <- Y
Yc[Y <= yc] <- NA
upper <- Y
upper[Y <= yc] <- yc 
lower <- Y
lower[Y <= yc] <- -10.0 
```

# Simulation naive

On crée un nouveau champ *Y.rn* de données randomisées: on conserve inchangées les valeurs supérieures à la coupure et on randomise les valeurs inférieures à la coupure.

```{r Randomisation}
pp <- sum(Y <= yc) / length(Y)
nc <- sum(Y <= yc)
ir <- sample(x = 1:nc, replace = FALSE, size = nc)
Y.rn <- Y
Y.rn[Y <= yc] <- qnorm(p = (ir)/nc * pp)
plot(x = Y, y = Y.rn, pch = 19, col = "lightblue", 
     xlim = 3.5 * c(-1,1), ylim = 3.5 * c(-1,1),
     main="Données randomisées vs. Données Initiales")
abline(a = 0, b = 1, col = "orange")
```

Les histogrammes de *Y* et *Y.rn*  sont bien identiques. 

```{r Histogrammes}
hist(Y, col = "gray",main="Histogramme des données initiales")
hist(Y.rn, col = "lightblue",main="Histogramme des données randomisées")
```

La base de données *data* est complétée. Elle comprend désormais:

- *Y*: la variable gaussienne initiale
- *Y.rn*: La variable gaussienne où les valeurs inférieures au seuil ont été randomisées

```{r Completion_Base_de_donnees}
data <- db.add(data, Y.rn, upper, lower)
data <- db.delete(data, names = "Y")
data <- db.rename(data, names = "Y.Y", newnames = "Y")
```

# Simulation avec Gibbs

Dans cette deuxième partie, les valeurs en-dessous du seuil sont simulées à l'aide d'un échantillonneur de Gibbs qui conserve la structure spatiale, tout en reconstituant des valeurs gaussiennes sous le seuil.

```{r Gibbs, echo = TRUE, eval = TRUE}
data <- db.locerase(data, loctype = "z")
data <- db.locate(data, names = "upper", loctype = "upper")
data <- db.locate(data, names = "lower", loctype = "lower")
data <- gibbs(db = data, model = m, seed = 232132, nbsimu = 1,
              nboot = 10, niter = 1000, flag.norm=FALSE,
              flag.ce = FALSE, radix = "Gibbs", modify.target = TRUE,
              verbose = FALSE)
```

Le nuage de corrélation entre la variable gaussienne initiale *Y* et la variable créée par Gibbs montre la validité de la procédure.

```{r Correlation}
dum = correlation(data, name1 = "Y", name2 = "Gibbs.G1", 
            flag.aspoint = TRUE, pch = 19, col = "lightblue", 
            flag.iso = TRUE, flag.diag = TRUE,
            title="Corrélation entre données Gibbs et données randomisées")
abline(h = yc)
```

# Comparaison

Il est maintenant intéressant de comparer les variogrammes expérimentaux de:

- *Y* la variable gaussienne initiale
- *Y.rn*  où les valeurs (sous le seuil) ont été obtenues par randomisation, sans respect de la structure spatiale
- *Gibbs.G1*: Où les valeurs sous le seuil ont été obtenues par Gibbs qui respecte la structure spatiale

```{r Variogrammes, echo = TRUE, eval = TRUE}
var_Y <- vario.calc(
  db.locate(data, names = "Y", loctype = "z", flag.locnew = TRUE),
  lag = 0.05, nlag = 10
  )

var_Y.rn <- vario.calc(
  db.locate(data, names = "Y.rn", loctype = "z", flag.locnew = TRUE),
  lag = 0.05, nlag = 10
  )

var_Y.gibbs <- vario.calc(
  db.locate(data, names = "Gibbs.G1", loctype = "z", flag.locnew = TRUE),
  lag = 0.05, nlag = 10
  )
```

```{r Variogrammes_individuels}
plot(var_Y, title = "Variable initiale")
plot(m, col = "red", add = TRUE)

plot(var_Y.rn, title = "Variable complétée par Randomisation")
plot(m, col = "red", add = TRUE)

plot(var_Y.gibbs, title = "Variable complétée par Gibbs")
plot(m, col = "red", add = TRUE)
```


```{r Superposition_variogrammes}
plot(var_Y, col="black", title = "Variable initiale (noir), randomisées (rouge) et Gibbs (blue)")
plot(m, col = "black", lwd=2, add = TRUE)
plot(var_Y.rn, col="red",add=TRUE)
plot(var_Y.gibbs, col="blue", add=TRUE)
```