Librerías y lectura de datos

Cargamos las librerías que utilizaremos durante la construcción del modelo, especificando las funciones que usamos en cada una de ellas.

library(lpSolve) # lp
library(readxl) # read_xlsx
library(ggplot2) # ggplot
library(caret) # createDataPartition

# Para que no aparezcan los resultados en notación científica
options("scipen"=100, digits=3)

Leemos los datos y visualizamos parte de ellos.

data_paises <- read_xlsx('data.xlsx', col_names = TRUE)
head(data_paises)
## # A tibble: 6 x 5
##   paises_oecd num_medicos num_camas gasto_salud_per_cap esp_vida
##   <chr>             <dbl>     <dbl>               <dbl>    <dbl>
## 1 Alemania           3.92       8.2               4750.     80.5
## 2 Austria            4.84       7.7               4966.     80.9
## 3 Bélgica            2.94       6.3               4543.     80.4
## 4 Canadá             2.40       2.7               5353.     81.6
## 5 Colombia           1.74       1.5                448.     75.9
## 6 Corea              2.08      10.3               1576.     80.8

Eliminamos la columna \(\textit{paises_oecd}\) porque no nos aporta información relevante para el presente análisis.

data <- data.frame(subset(data_paises, select = -paises_oecd))

Definición de las variables y los parámetros

Definimos las variables conocidas, el parámetro \(\varepsilon\) y las dimensiones \(d\) y \(n\). Además, se dividen los datos en entrenamiento y testeo.

# División datos entrenamiento/testeo
inTrain <- createDataPartition(y=data$esp_vida, p=0.8, list=FALSE)
train <- data[inTrain,]
test <- data[-inTrain,]

# Variables conocidas
x <- train[, 1:ncol(train)-1]
y <- train[, ncol(train)]

d <- ncol(x) # Dimensión de las variables x
n <- nrow(x) # Tamaño de la muestra

epsilon <- 0.5 # Margen de error permitido

Construcción del problema de optimización (SVR2)

A continuación, vamos a definir el vector correspondiente a la función objetivo \(obj\), y la matriz de coeficientes \(A\) y el vector de direcciones \(dir\) correspondientes a las restricciones.

Es importante tener en cuenta que el orden de las variables elegido es el siguiente:

\[ x = (w_1^{+} \quad ... \quad w_d^{+} \quad w_1^{-} \quad ... \quad w_d^{-} \quad b^{+} \quad b^{-} \quad \xi_1 \quad ... \quad \xi_n \quad \xi_1^{\ast} \quad ... \quad \xi_n^{\ast}) \]

Podemos ver que a la variable \(b \in \mathbb{R}\) también le hemos aplicado la transformación del tipo \(b:=b^{+}-b^{-}\), con \(b^{+}, b^{-} \geq 0\). Esto es debido a que la función \(lp\) del paquete ‘lpSolve’ da por hecho que todas las variables son no negativas.

obj <- c(rep(1,2*d), rep(0,2+2*n))

Restricción (1): \(\sum\nolimits_{i=1}^{n}\left( \xi _{i}+\xi _{i}^{\ast }\right) \leq c\)

A1 <- c(rep(0,2*d+2),rep(1,2*n))

Restricción (2): \(y_{i}-\left( \left\langle \omega ^{+}-\omega ^{-},x_{i}\right\rangle + (b^{+} - b^{-})\right) \leq \varepsilon +\xi _{i},~~i=1,...,n\)

A2 <- cbind(-x,x,matrix(-1,n,1),matrix(1,n,1),-diag(n),matrix(0,n,n))

Restricción (3): \(\left( \left\langle \omega ^{+}-\omega ^{-},x_{i}\right\rangle + (b^{+} - b^{-})\right) - y_{i} \leq \varepsilon +\xi _{i},~~i=1,...,n\)

A3 <- cbind(x,-x,matrix(1,n,1),matrix(-1,n,1),matrix(0,n,n),-diag(n))

Unimos las restricciones \(A_i, i = 1,2,3\) para formar la matriz final A:

colnames(A2) <- 1:ncol(A2)
colnames(A3) <- 1:ncol(A2)

A <- rbind(A1,A2,A3)
dir <- rep("<=", 1+2*n)

Nos falta definir el vector \(b\), en el cual interviene el parámetro \(c>0\), que es la cota de la suma de los errores (ver Restricción (1)).

Para ello, se construye el siguiente bucle, donde vamos a ir dándole valores al parámetro \(c\) desde 1 hasta un valor grande y guardando las soluciones VÁLIDAS en la matriz \(solutionC\), en función de dos condiciones:

  1. Condición para que comience a almacenar las soluciones en la matriz \(solutionC\): el problema de optimización sea FACTIBLE. Para ello, se tiene que dar que el estado de la solución (\(sol\$status\)) tome el valor 0. Por lo tanto, el bucle comenzará dándole valores a \(c = 1, 2, 3, ...\), para los cuales el estado de la solución sea no factible, hasta que llegue a un valor de \(c\) para el cual sí sea factible, y entonces se comenzará a guardar dicha solución.

  2. Condición de parada: |sol$objval - last$objval| < 0.000001. La diferencia, en valor absoluto, del valor de la función objetivo en la iteración (i) y el valor de la función objetivo en la iteración anterior (i-1) ha de ser lo suficientemente pequeña (0.000001) como para que ya no sea de interés seguir almacenando las soluciones para los diferentes valores de \(c\). En el bucle lo definimos al revés, es decir, mientras la diferencia sea mayor que ese valor muy pequeño, que guarde la solución.

solutionC <- vector() # Vector vacío para almacenar w+, w- y b
chis <- vector() # Vector vacío para almacenar chi, chi*
last_objval <- .Machine$integer.max # Comparativa primera iteración, un valor muy grande cualquiera (2147483647)

for (i in 1:last_objval){ # Valores para c
  b <- matrix(cbind(c(i, epsilon-y, epsilon+y)), ncol=1) # Vector términos independientes b
  sol <- lp("min", obj, A, dir, b) # Resuelvo y guardo en 'sol' la solución
  if (sol$status == 0){ # Primera condición
    if (abs(sol$objval-last_objval)>0.000001) # Segunda condición a la inversa
    {
      solutionC <- rbind(solutionC, c(i, sol$solution[1:(2*d+2)])) # Almaceno el valor de 'c' en la primera columna y la solución de w+,w- y b
      last_objval <- sol$objval # Guardo en 'last_objval' el valor de la f. obj de la iteración actual
      chis <- rbind(chis, c(sol$solution[(2*d+2+1):length(sol$solution)])) # Almaceno en 'chis' los valores de chi,chi*
    } else {break} # Parada cuando deje de cumplirse la Segunda condición a la inversa
  }
}

# head(solutionC) # Solución variables w+, w- y b
# head(chis) # Solución variables chi, chi*

De esta forma, ya tenemos las soluciones válidas para los diferentes valores de \(c\) almacenadas en el vector \(solutionC\). Ahora, procedemos a deshacer las transformaciones realizadas anteriormente a las variables \(\omega=\omega^{+}-\omega^{-}\), \(b=b^{+}-b^{-}\), y guardar las soluciones definitivas de \(w,b\) en la matriz \(solution\_fin\), y \(\xi,\xi^{*}\) en \(chi\) y \(chi\_ast\), respectivamente.

# Almaceno en solution_fin los valores de w, b con su correspondiente valor c
solution_fin <- vector() # Vector vacío para w, b

for (i in 1:nrow(solutionC)){
  vector <- solutionC[i,1] # Columna valor c
  for (j in 2:(d+1)){
    vector <- cbind(vector, solutionC[i,j]-solutionC[i,j+d]) # Columnas valores w+ - w-
  }
  vector <- cbind(vector, solutionC[i,ncol(solutionC)-1]-solutionC[i,ncol(solutionC)]) # Columna valor b
  solution_fin <- rbind(solution_fin,vector)
}
solution_fin <- data.frame(matrix(solution_fin,ncol=1+d+1))
names(solution_fin) <- c("c","w1","w2","w3","b")
head(solution_fin)
##    c    w1 w2       w3    b
## 1 38 0.186  0 0.000728 76.9
## 2 39 0.000  0 0.000645 77.9
## 3 40 0.000  0 0.000423 78.9
## 4 41 0.000  0 0.000303 79.3
## 5 42 0.000  0 0.000234 79.7
## 6 43 0.000  0 0.000183 80.0
# Almaceno en chi y chi_ast sus correspondientes valores
chi <- chis[,1:(n)]
chi_ast <- chis[,(n+1):(ncol(chis))]

Gráfico \(\sum\nolimits_{i=1}^{d}\left\vert \omega _{i}\right\vert\) frente valores \(c\)

En este gráfico vamos a ver cómo evoluciona la suma de los valores de \(\omega\) en valor absoluto, en función de los valores de \(c\).

# Guardo en data.frames los correspondientes valores de la suma(w), b y el parámetro c
w <- solution_fin[,2:(ncol(solution_fin)-1)]
sum_w <- data.frame(sum_w=rowSums(abs(w)))
b <- data.frame(b=solution_fin[,ncol(solution_fin)])
c <- data.frame(c=solution_fin[,1])
sum_w_c <- cbind(sum_w,c)

# Gráfico
ggplot(sum_w_c, aes(x=c, y=sum_w)) +
  geom_point() +
  labs(title="Suma valores absolutos de w en función de c", x="Valor c", y="Suma 'abs(w)'")

A mayores valores de \(c\), menor es la suma en valor absoluto de \(\omega\). En otras palabras, cuanto menos restrictivos somos con la restricción (1), más horizontal será el tubo que envuelve a los puntos.

Validación del modelo

Sobre los datos de entrenamiento

Para ver si el modelo es coherente con los datos que se han utilizado para “entrenarlo”, vamos a comparar la diferencia entre real y estimado con los valores de \(\xi,\xi^{\ast}\) sumándoles \(\varepsilon\). Es decir;

\[ \text{Comparar } |\hat{y}-y| \text{ frente a } \xi+\varepsilon,\xi^{\ast}+\varepsilon \] La comparativa la vamos a hacer para dos casos: la solución válida para el valor de \(c\) más pequeño y para el más grande.

Caso 1: hiperplano \(c\) más pequeño

Guardo la solución para el valor de \(c\) más pequeño.

w_C_min <- as.matrix(w[1,])
b_min <- as.matrix(b[1,])
chi_min <- as.matrix(chi[1,])
chi_ast_min <- as.matrix(chi_ast[1,])
chi_chi_ast <- cbind(chi_min,chi_ast_min)
sum_chis <- chi_min+chi_ast_min

Calculamos \(\hat{y} = \left\langle \hat{\omega} ,x_{i}\right\rangle + \hat{b}, \quad i=1,...,n\), para los datos de entrenamiento.

y_est_1 <- vector()
for (i in 1:nrow(x)){
  y_est_1 <- rbind(y_est_1, (w_C_min%*%t(x)[,i]) + b_min)
}

A continuación, calculamos la diferencia entre estimado y real, \(|\hat{y}-y|\).

dif_est_real_1 <- abs(y_est_1-as.matrix(y))

Por último, hacemos la comparación \(|\hat{y}-y| \text{ frente a } \xi+\varepsilon,\xi^{\ast}+\varepsilon\).

comparativa_1 <- data.frame(cbind(sum_chis+epsilon, dif_est_real_1), row.names=NULL)
names(comparativa_1) <- c("Chi_Mas_Epsilon","Est_Menos_Real"); comparativa_1
##    Chi_Mas_Epsilon Est_Menos_Real
## 1            0.534         0.5341
## 2            0.500         0.4664
## 3            0.500         0.4193
## 4            2.397         2.3973
## 5            1.327         1.3274
## 6            2.945         2.9447
## 7            4.754         4.7540
## 8            1.908         1.9077
## 9            0.500         0.0387
## 10           1.106         1.1062
## 11           1.183         1.1834
## 12           3.101         3.1009
## 13           0.500         0.3901
## 14           2.629         2.6289
## 15           2.459         2.4587
## 16           2.359         2.3587
## 17           1.992         1.9918
## 18           4.228         4.2276
## 19           4.448         4.4483
## 20           2.713         2.7126
## 21           2.753         2.7533
## 22           1.147         1.1472
## 23           1.332         1.3316
## 24           0.965         0.9652
## 25           0.500         0.5000
## 26           0.500         0.5000
## 27           2.220         2.2204
dif_comparativa_1 <- as.data.frame(comparativa_1[,1]-comparativa_1[,2])
names(dif_comparativa_1) <- c("Diferencias"); round(dif_comparativa_1, digits=3)
##    Diferencias
## 1        0.000
## 2        0.034
## 3        0.081
## 4        0.000
## 5        0.000
## 6        0.000
## 7        0.000
## 8        0.000
## 9        0.461
## 10       0.000
## 11       0.000
## 12       0.000
## 13       0.110
## 14       0.000
## 15       0.000
## 16       0.000
## 17       0.000
## 18       0.000
## 19       0.000
## 20       0.000
## 21       0.000
## 22       0.000
## 23       0.000
## 24       0.000
## 25       0.000
## 26       0.000
## 27       0.000

Se puede observar que se cumplen las restricciones que se han definido en el modelo.

Caso 2: hiperplano \(c\) más grande

En este punto hacemos lo mismo que en el anterior, pero para el valor de \(c\) máximo dentro de nuestra matriz de soluciones.

w_C_max <- as.matrix(w[nrow(w),])
b_max <- as.matrix(b[nrow(b),])
chi_max <- as.matrix(chi[nrow(w),])
chi_ast_max <- as.matrix(chi_ast[nrow(w),])
chi_chi_ast <- cbind(chi_max,chi_ast_max)
sum_chis <- chi_max+chi_ast_max
y_est_2 <- vector()
for (i in 1:nrow(x)){
  y_est_2 <- rbind(y_est_2, (w_C_max%*%t(x)[,i]) + b_max)
}
dif_est_real_2 <- abs(y_est_2-as.matrix(y))
comparativa_2 <- data.frame(cbind(sum_chis+epsilon, dif_est_real_2), row.names=NULL)
names(comparativa_2) <- c("Chi_Epsilon","Est_Real"); comparativa_2
##    Chi_Epsilon Est_Real
## 1        0.500   0.1927
## 2        0.590   0.5902
## 3        1.302   1.3024
## 4        0.500   0.4732
## 5        0.500   0.2220
## 6        2.080   2.0805
## 7        1.605   1.6049
## 8        4.020   4.0195
## 9        0.500   0.2805
## 10       1.622   1.6220
## 11       0.500   0.2878
## 12       5.283   5.2829
## 13       0.500   0.5000
## 14       2.571   2.5707
## 15       1.359   1.3585
## 16       1.893   1.8927
## 17       2.750   2.7498
## 18       6.568   6.5683
## 19       6.483   6.4829
## 20       5.380   5.3803
## 21       1.733   1.1049
## 22       3.600   3.6000
## 23       0.500   0.0268
## 24       0.559   0.5585
## 25       2.271   2.2707
## 26       1.359   1.3585
## 27       4.973   4.9733
dif_comparativa_2 <- as.data.frame(comparativa_2[,1]-comparativa_2[,2])
names(dif_comparativa_2) <- c("Diferencias"); round(dif_comparativa_2,digits=3)
##    Diferencias
## 1        0.307
## 2        0.000
## 3        0.000
## 4        0.027
## 5        0.278
## 6        0.000
## 7        0.000
## 8        0.000
## 9        0.220
## 10       0.000
## 11       0.212
## 12       0.000
## 13       0.000
## 14       0.000
## 15       0.000
## 16       0.000
## 17       0.000
## 18       0.000
## 19       0.000
## 20       0.000
## 21       0.628
## 22       0.000
## 23       0.473
## 24       0.000
## 25       0.000
## 26       0.000
## 27       0.000

De la misma forma, se cumplen las restricciones que se han definido en el modelo.

Sobre los datos de testeo

En esta última sección vamos a aplicar el modelo sobre los datos de testeo y a calcular una medida de bondad de ajuste para el mismo, definida en el TFG como \(MBA\).

\[ MBA = \left(\frac{1}{n}\sum\nolimits_{i=1}^{n} | y_{i}-\hat{y}_i |\right) \frac{1}{\frac{1}{n}\sum\nolimits_{i=1}^{n} | y_{i}-\overline{y} |} \]

Lo haremos, como en la sección anterior, para el caso de \(c\) más pequeño y para el más grande.

x_test <- test[,1:ncol(test)-1]
y_test <- test[,ncol(test)]

Caso 1: hiperplano \(c\) más pequeño

Calculamos \(\hat{y} = \left\langle \hat{\omega} ,x_{i}\right\rangle + \hat{b}, \quad i=1,...,n\), para los datos de testeo.

y_est_1 <- vector()
for (i in 1:nrow(x_test)){
  y_est_1 <- rbind(y_est_1, w_C_min%*%t(x_test)[,i] + b_min)
}

Y aplicamos la Medida de Bondad de Ajuste del modelo.

MBA_1 <- ((mean(abs(y_test-y_est_1)))/(mean(abs(y_test-mean(y_test))))); MBA_1
## [1] 0.578

Caso 2: hiperplano \(c\) más grande

Realizamos los mismos pasos que en el punto anterior.

y_est_2 <- vector()
for (i in 1:nrow(x_test)){
  y_est_2 <- rbind(y_est_2, w_C_max%*%t(x_test)[,i] + b_max)
}
MBA_2 <- ((mean(abs(y_test-y_est_2)))/(mean(abs(y_test-mean(y_test))))); MBA_2
## [1] 0.939