Modelos de regresión ambiental

Análisis espacial con R

Dr. Francisco Zambrano

Preparar datos para interpolación

library(tidyverse)
library(sf)
library(tmap)
library(terra)

region <- read_rds('data/regs_utm_sf.rds')
data_temp <- read_rds('data/data_temp_enero.rds')
preds_ene <- rast('data/predictores_enero.tif')
grilla <- preds_ene[[1]]
values(grilla) <- NA
set.seed(234)
data_temp_mod <- data_temp |> sample_frac(0.8)
data_temp_val <- data_temp |> filter(!(ema %in% data_temp_mod$ema))

calc_met <- function(ras,object_sf,var = 'temp_enero',k=1) {
  df_val <- terra::extract(ras,terra::vect(object_sf)) |> 
    cbind(object_sf[,var]) 
  rmse <- function(obs,pred) sqrt(sum((obs-pred)^2)/length(obs))
  rmse_v <- rmse(df_val[,3],df_val[,2])
  r2 <- cor(df_val[,3],df_val[,2])^2
  r2_adj <- 1- (1-r2)*(dim(df_val)[1]-1)/(dim(df_val)[1]-k-1) 
  c('rmse' = rmse_v,'rsq' = r2,'rsq_adj' = r2_adj)
}

Metricas para evaluación

Raíz cuadrada del error cuadratico medio

\[RMSE = \sqrt{\frac{\sum{(obs-est)^2}}{N}}\] Coeficiente de determinación R cuadrado

\[R^2=1-\frac{\sum{(y-\hat{y})^2}}{\sum(y-\bar{y})^2}=1-\frac{Variación\,no\,explicada}{Variación\,total}\]

Regresión lineal

\[\hat Z(s_0) = \hat \beta_0 + \beta_1\cdot X_1(s_0) + ...+\hat \beta_p\cdot q_p(s_0) = \sum_{k=0}^p \hat \beta_k \cdot X_k(s_0)\]

Podemos ajustar el modelo de regresión utilizando la función lm (regresion simple), glm (regresión generalizada). Luego utilizar terra::predict para predecir en el raster.

mod_lm1 <- lm(temp_enero~dem,data_temp_mod)
temp_lm1 <- terra::predict(preds_ene,mod_lm1,se.fit = TRUE)
(met_lm1 <- calc_met(temp_lm1,data_temp_val))

      rmse        rsq    rsq_adj 
18.9753110  0.9333186  0.9320362

Regresión lineal

Utilizando la función krige de {gstat}, utiliza modelos de regresión lineal generalizada en vez de modelos de regresión simple.

library(stars)
library(gstat)
preds_ene_st <- preds_ene |> st_as_stars() |>  split('band')
temp_lm2 <- krige(temp_enero~dem,data_temp_mod,preds_ene_st)

[ordinary or weighted least squares prediction]

temp_lm2 <- rast(temp_lm2)
(met_lm2 <- calc_met(temp_lm2[[1]],data_temp_val))

     rmse       rsq   rsq_adj 
3.9278719 0.4871318 0.4772689

Regresión lineal

Si incorporamos más predictores

temp_lm3 <- krige(temp_enero~dem+lst_ene,data_temp_mod,preds_ene_st)

[ordinary or weighted least squares prediction]

temp_lm3 <- rast(temp_lm3)
(met_lm3 <- calc_met(temp_lm3[[1]],data_temp_val,k=2))

     rmse       rsq   rsq_adj 
2.8348459 0.6142652 0.5991384

Regresión lineal

… y más predictores

temp_lm4 <- krige(temp_enero~dem+lst_ene+dist_costa+ndvi_ene,data_temp_mod,preds_ene_st)

[ordinary or weighted least squares prediction]

temp_lm4 <- rast(temp_lm4)
(met_lm4 <- calc_met(temp_lm4[[1]],data_temp_val,k=4))

     rmse       rsq   rsq_adj 
2.4003937 0.6142051 0.5827117

Regresión lineal

Ahora uitilizamos los seis predictores.

temp_lm5 <- krige(temp_enero~dem+aspect+slope+dist_costa+ndvi_ene+lst_ene,data_temp_mod,preds_ene_st)

[ordinary or weighted least squares prediction]

temp_lm5 <- temp_lm5 |> rast()
(met_lm5 <- calc_met(temp_lm5[[1]],data_temp_val,k=6))

     rmse       rsq   rsq_adj 
2.3682991 0.6231314 0.5750205

Evaluación de los modelos

Regresión lineal

Varianza del error

Hay una diferencia entre los intervalos de confianza respecto a la linea de regresión en el \(x_0\) y la de predicción de \(y\) en el punto \(x_0\).
Linea \(x_0\) se,

\[\hat \sigma\sqrt{\frac{1}{n} + \frac{(x_0 - \bar X)^2}{\sum_{i=1}^n (X_i - \bar X)^2}}\]

Intervalos de predicción se en \(x_0\),

\[\hat \sigma\sqrt{1 + \frac{1}{n} + \frac{(x_0 - \bar X)^2}{\sum_{i=1}^n (X_i - \bar X)^2}}\]

Regresion Lineal

Evaluación del modelo: Validación cruzada

Leave-one-out

Regresion Lineal

Evaluación del modelo: Validación cruzada

K-Folds

Regresion Lineal

Evaluación del modelo: Validación cruzada

Leave-one-out

(temp_lm5_loocv <- krige.cv(temp_enero~dem+aspect+slope+dist_costa+ndvi_ene+lst_ene,data_temp_mod))

Simple feature collection with 218 features and 6 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 106716.7 ymin: 5633261 xmax: 400167.8 ymax: 6748953
Projected CRS: WGS 84 / UTM zone 19S
First 10 features:
   var1.pred var1.var observed     residual       zscore fold
1   21.47829 2.927101 21.19634 -0.281957714 -0.164803011    1
2   20.11498 2.861589 20.46371  0.348733893  0.206153419    2
3   20.33654 2.862717 20.33865  0.002113285  0.001249018    3
4   19.52156 2.850767 17.84503 -1.676535973 -0.992960512    4
5   20.83129 2.956989 20.38382 -0.447470941 -0.260219598    5
6   22.05254 3.173598 21.71897 -0.333570747 -0.187245751    6
7   20.00533 2.841105 21.49471  1.489377321  0.883610990    7
8   18.54163 2.887854 18.21565 -0.325983188 -0.191826048    8
9   20.04311 2.849239 20.26546  0.222342898  0.131722151    9
10  19.50606 2.892990 18.53133 -0.974725849 -0.573071790   10
                   geometry
1  POINT (316866.4 6196602)
2  POINT (342319.3 6268991)
3    POINT (342450 6269810)
4  POINT (310362.3 6505139)
5  POINT (260348.8 6189793)
6  POINT (331096.5 6180897)
7  POINT (274890.3 6090475)
8    POINT (229287 6111926)
9  POINT (325277.7 6262424)
10 POINT (249767.1 6210853)

Regresion Lineal

Evaluación del modelo: Validación cruzada

K-Fold

(temp_lm5_kfoldcv <- krige.cv(temp_enero~dem+aspect+slope+dist_costa+ndvi_ene+lst_ene,data_temp_mod,nfold=10))

Simple feature collection with 218 features and 6 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 106716.7 ymin: 5633261 xmax: 400167.8 ymax: 6748953
Projected CRS: WGS 84 / UTM zone 19S
First 10 features:
   var1.pred var1.var observed    residual      zscore fold
1   21.61702 2.971277 21.19634 -0.42067867 -0.24405008    9
2   20.12711 3.026410 20.46371  0.33659919  0.19348583   10
3   20.31389 2.700538 20.33865  0.02476331  0.01506897    8
4   19.52877 3.027277 17.84503 -1.68374342 -0.96772029   10
5   20.64375 2.994194 20.38382 -0.25992841 -0.15021516    7
6   21.95137 3.236042 21.71897 -0.23240105 -0.12919062    7
7   20.00204 3.011960 21.49471  1.49266788  0.86007945    2
8   18.36952 2.883719 18.21565 -0.15387331 -0.09061223    3
9   20.05023 3.008518 20.26546  0.21522407  0.12408365    2
10  19.43614 2.918717 18.53133 -0.90480693 -0.52961458    7
                   geometry
1  POINT (316866.4 6196602)
2  POINT (342319.3 6268991)
3    POINT (342450 6269810)
4  POINT (310362.3 6505139)
5  POINT (260348.8 6189793)
6  POINT (331096.5 6180897)
7  POINT (274890.3 6090475)
8    POINT (229287 6111926)
9  POINT (325277.7 6262424)
10 POINT (249767.1 6210853)