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(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.8Eliminamos 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))Definimos las variables conocidas, el parámetro \(\varepsilon\) y las dimensiones \(d\) y \(n\). Además, vamos a dividir 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 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(0,2*d+2),rep(1,2*n))Restricción (1): \(\sum\nolimits_{i=1}^{d}\left( \omega_{i}^{+}+\omega_{i}^{-}\right)\leq h\)
A1 <- c(rep(1,2*d),rep(0,2+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 \(h>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 \(h\) desde 0.1 hasta un valor grande y guardando las soluciones VÁLIDAS en la matriz \(solutionH\), en función de la siguiente condición:
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.
solutionH <- vector() # Vector vacío para almacenar w+, w- y b
chis <- vector() # Vector vacío para almacenar chi, chi*
last_objval <- 10000 # Comparativa primera iteración frente a un valor muy grande
for (i in seq(0.1,last_objval,0.1)){ # Valores para h
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 (abs(sol$objval-last_objval)>0.000001) # Condición de parada a la inversa
{
solutionH <- rbind(solutionH,c(i,sol$solution[1:(2*d+2)])) # Almaceno el valor de 'h' 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 condición inversa
}
head(solutionH) # Solución variables w+, w- y b## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [1,] 0.1 0.000 0.0994 0.000622 0 0 0 77.0 0
## [2,] 0.2 0.000 0.1994 0.000644 0 0 0 76.3 0
## [3,] 0.3 0.127 0.1721 0.000646 0 0 0 76.1 0
## [4,] 0.4 0.231 0.1679 0.000660 0 0 0 75.7 0
## [5,] 0.5 0.344 0.1554 0.000641 0 0 0 75.4 0
## [6,] 0.6 0.450 0.1496 0.000622 0 0 0 75.1 0# head(chis) # Solución variables chi, chi*De esta forma, ya tenemos las soluciones válidas para los diferentes valores de \(h\) almacenadas en el vector \(solutionH\). 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 h
solution_fin <- vector() # Vector vacío para w, b
for (i in 1:nrow(solutionH)){
vector <- solutionH[i,1] # Columna valor h
for (j in 2:(d+1)){
vector <- cbind(vector, solutionH[i,j]-solutionH[i,j+d]) # Columnas valores w+ - w-
}
vector <- cbind(vector, solutionH[i,ncol(solutionH)-1]-solutionH[i,ncol(solutionH)]) # Columna valor b
solution_fin <- rbind(solution_fin,vector)
}
solution_fin <- data.frame(matrix(solution_fin,ncol=1+d+1))
names(solution_fin) <- c("h","w1","w2","w3","b")
head(solution_fin)## h w1 w2 w3 b
## 1 0.1 0.000 0.0994 0.000622 77.0
## 2 0.2 0.000 0.1994 0.000644 76.3
## 3 0.3 0.127 0.1721 0.000646 76.1
## 4 0.4 0.231 0.1679 0.000660 75.7
## 5 0.5 0.344 0.1554 0.000641 75.4
## 6 0.6 0.450 0.1496 0.000622 75.1# Almaceno en chi y chi_ast sus correspondientes valores
chi <- chis[,1:(n)]
chi_ast <- chis[,(n+1):(ncol(chis))]
# Guardo en diferentes vectores las variables calculadas y el parámetro h
w <- solution_fin[,2:(ncol(solution_fin)-1)]
b <- data.frame(b=solution_fin[,ncol(solution_fin)])
h <- data.frame(c=solution_fin[,1])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 \(h\) más pequeño y para el más grande.
Guardo la solución para el valor de \(h\) 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_minCalculamos \(\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.500 0.2758
## 2 0.500 0.0373
## 3 0.500 0.1117
## 4 1.006 1.0062
## 5 1.593 1.5931
## 6 1.769 1.7691
## 7 1.402 1.4022
## 8 3.472 3.4716
## 9 3.839 3.8387
## 10 1.897 1.8970
## 11 1.385 1.3846
## 12 1.926 1.9260
## 13 3.277 3.2771
## 14 2.866 2.8659
## 15 2.911 2.9111
## 16 1.478 1.4776
## 17 4.319 4.3191
## 18 4.477 4.4769
## 19 0.582 0.5822
## 20 2.580 2.5804
## 21 1.569 1.5687
## 22 1.795 1.7954
## 23 1.406 1.4062
## 24 0.500 0.5000
## 25 0.500 0.5000
## 26 0.599 0.5986
## 27 2.269 2.2690dif_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.224
## 2 0.463
## 3 0.388
## 4 0.000
## 5 0.000
## 6 0.000
## 7 0.000
## 8 0.000
## 9 0.000
## 10 0.000
## 11 0.000
## 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.000
## 22 0.000
## 23 0.000
## 24 0.000
## 25 0.000
## 26 0.000
## 27 0.000Se puede observar que se cumplen las restricciones que se han definido en el modelo.
En este punto hacemos lo mismo que en el anterior, pero para el valor de \(h\) 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_maxy_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.5000
## 2 0.688 0.6881
## 3 0.500 0.2827
## 4 1.880 1.8796
## 5 0.500 0.5000
## 6 2.378 2.3777
## 7 1.990 1.9897
## 8 3.413 3.4134
## 9 2.909 2.9090
## 10 1.808 1.8079
## 11 1.644 1.6443
## 12 0.500 0.5000
## 13 3.148 3.1484
## 14 3.065 3.0655
## 15 2.834 2.8342
## 16 2.008 2.0081
## 17 4.152 4.1521
## 18 4.961 4.9606
## 19 0.500 0.0509
## 20 1.659 1.6588
## 21 1.645 1.6453
## 22 1.548 1.5476
## 23 2.012 2.0117
## 24 0.652 0.6519
## 25 0.500 0.4797
## 26 0.500 0.5000
## 27 1.219 1.2185dif_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.000
## 2 0.000
## 3 0.217
## 4 0.000
## 5 0.000
## 6 0.000
## 7 0.000
## 8 0.000
## 9 0.000
## 10 0.000
## 11 0.000
## 12 0.000
## 13 0.000
## 14 0.000
## 15 0.000
## 16 0.000
## 17 0.000
## 18 0.000
## 19 0.449
## 20 0.000
## 21 0.000
## 22 0.000
## 23 0.000
## 24 0.000
## 25 0.020
## 26 0.000
## 27 0.000De la misma forma, se cumplen las restricciones que se han definido en el modelo.
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 \(h\) más pequeño y para el más grande.
x_test <- test[,1:ncol(test)-1]
y_test <- test[,ncol(test)]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.734Realizamos 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.816