El objetivo del presente documento es el de crear un breve informe con RStudio mediante el análisis numérico de un conjunto de datos, en este caso, regresión múltiple con variable dependiente discreta.

Tras abrir RStudio vamos a crear una carpeta y establecerla como el directorio de trabajo (Session -> Set Working Directory -> Choose Directoy), abrimos un nuevo proyecto de R Markdown (muy interesante consultar su ayuda: Help -> Markdown Quick Reference). Seleccionaremos que genere un html.

Importación de datos

El primer paso será incorporar la información a tratar. Puesto que normalmente se disponen los datos en una hoja excel, los guardaremos como un archivo .csv separado por puntos y comas (como es este caso) o como un archivo .txt separado por tabulaciones. Los datos se pueden incorporar a RStudio siguiendo el camino Tools -> Import Dataset o mediante el siguiente comando (se supone que los datos se encuentran en la carpeta que hemos establecido como directorio de trabajo):

datos = read.table("credito.csv", header=T, sep=";")
head(datos) # vemos los datos, en la primera línea están el nombre de las variables
##   credito ingresos laboral cargas
## 1       1      4.0       2      0
## 2       1      4.5       2      1
## 3       1      5.0       2      0
## 4       1      4.0       2      0
## 5       0      2.5       1      1
## 6       0      2.5       1      1
# credito   = devolución del crédito (1 si el cliente devolvió el crédito)                  
# ingresos = ingresos anuales brutos del cliente (en decenas de miles de euros)                     
# laboral   = situación labotal del cliente (0 en paro, 1 contrato temporal, 2 contrato fijo)                       
# cargas = cargas familiares (1 con cargas familiares)                      
attach(datos) # incorporamos los datos a R para referenciar cada variable por el nombre anterior

Algunas características de los datos son las siguientes:

table(credito,laboral)
##        laboral
## credito  0  1  2
##       0 15 12  1
##       1  0  9 20
table(credito,cargas)
##        cargas
## credito  0  1
##       0  9 19
##       1 22  7
summary(ingresos)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.40    2.50    4.00    6.17    8.20   25.00
tapply(ingresos,as.factor(credito),summary)
## $`0`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.400   2.000   2.500   3.118   3.700  12.000 
## 
## $`1`
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.000   6.000   8.000   9.117  12.000  25.000

Se observa, por ejemplo, que de los 29 individuos que devuelven el crédito, 22 no tienen cargas familiares, 20 tienen contrato indefinido y tienen unos ingresos medios de 91170 euros. Por otro lado, de los 28 que no lo devuelven, 19 tienen cargas familiares, sólo 1 tiene contrato indefinido y los ingresos anuales medios son de 31180 euros.

Si los datos se encuentran en una vez se podrían importar directamente mediante el código:

datosWeb <- read.table("http://www.ugr.es/local/romansg/material/WebEco/04-Eco2/Ordenador/R/credito.csv", header=T, sep=";")
head(datosWeb)  
##   credito ingresos laboral cargas
## 1       1      4.0       2      0
## 2       1      4.5       2      1
## 3       1      5.0       2      0
## 4       1      4.0       2      0
## 5       0      2.5       1      1
## 6       0      2.5       1      1

Regresión lineal múltiple

Si estamos interesados en analizar la influencia que tiene sobre la devolución de un crédito el resto de variables, por ejemplo, podemos realizar una regresión lineal múltiple estimada por Mínimos Cuadrados Ordinarios (MCO). Para ello recurrimos al código lm:

mlp = lm(credito~ingresos+as.factor(laboral)+cargas)
summary(mlp)
## 
## Call:
## lm(formula = credito ~ ingresos + as.factor(laboral) + cargas)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.74378 -0.14089  0.01514  0.13246  0.55486 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -0.006364   0.086622  -0.073   0.9417    
## ingresos             0.045123   0.008943   5.046 5.87e-06 ***
## as.factor(laboral)1  0.203331   0.099414   2.045   0.0459 *  
## as.factor(laboral)2  0.648285   0.103144   6.285 6.80e-08 ***
## cargas              -0.168883   0.075703  -2.231   0.0300 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2646 on 52 degrees of freedom
## Multiple R-squared:  0.7444, Adjusted R-squared:  0.7248 
## F-statistic: 37.87 on 4 and 52 DF,  p-value: 8.003e-15

Es interesante advertir cómo se han creado variables binarias a partir de la variable referente a la situación laboral del cliente y que la categoría de referencia es aquella en la que el cliente está en paro.

Otra opción es la de generar las variables dummys:

n = length(ingresos)
laboral0 = array(0,n)
laboral1 = array(0,n)
laboral2 = array(0,n)
for (i in 1:n){
  if(laboral[i] == 0){laboral0[i] = 1}
  if(laboral[i] == 1){laboral1[i] = 1}
  if(laboral[i] == 2){laboral2[i] = 1}
}
cbind(laboral, laboral0, laboral1, laboral2)
##       laboral laboral0 laboral1 laboral2
##  [1,]       2        0        0        1
##  [2,]       2        0        0        1
##  [3,]       2        0        0        1
##  [4,]       2        0        0        1
##  [5,]       1        0        1        0
##  [6,]       1        0        1        0
##  [7,]       2        0        0        1
##  [8,]       2        0        0        1
##  [9,]       0        1        0        0
## [10,]       0        1        0        0
## [11,]       2        0        0        1
## [12,]       1        0        1        0
## [13,]       2        0        0        1
## [14,]       2        0        0        1
## [15,]       0        1        0        0
## [16,]       0        1        0        0
## [17,]       0        1        0        0
## [18,]       1        0        1        0
## [19,]       2        0        0        1
## [20,]       1        0        1        0
## [21,]       1        0        1        0
## [22,]       1        0        1        0
## [23,]       0        1        0        0
## [24,]       1        0        1        0
## [25,]       0        1        0        0
## [26,]       1        0        1        0
## [27,]       1        0        1        0
## [28,]       2        0        0        1
## [29,]       1        0        1        0
## [30,]       1        0        1        0
## [31,]       0        1        0        0
## [32,]       2        0        0        1
## [33,]       1        0        1        0
## [34,]       2        0        0        1
## [35,]       2        0        0        1
## [36,]       1        0        1        0
## [37,]       2        0        0        1
## [38,]       1        0        1        0
## [39,]       0        1        0        0
## [40,]       1        0        1        0
## [41,]       1        0        1        0
## [42,]       2        0        0        1
## [43,]       1        0        1        0
## [44,]       2        0        0        1
## [45,]       0        1        0        0
## [46,]       0        1        0        0
## [47,]       1        0        1        0
## [48,]       2        0        0        1
## [49,]       0        1        0        0
## [50,]       1        0        1        0
## [51,]       2        0        0        1
## [52,]       0        1        0        0
## [53,]       0        1        0        0
## [54,]       0        1        0        0
## [55,]       2        0        0        1
## [56,]       2        0        0        1
## [57,]       1        0        1        0
mlp.bis = lm(credito~ingresos+laboral1+laboral2+cargas)
summary(mlp.bis)
## 
## Call:
## lm(formula = credito ~ ingresos + laboral1 + laboral2 + cargas)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.74378 -0.14089  0.01514  0.13246  0.55486 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.006364   0.086622  -0.073   0.9417    
## ingresos     0.045123   0.008943   5.046 5.87e-06 ***
## laboral1     0.203331   0.099414   2.045   0.0459 *  
## laboral2     0.648285   0.103144   6.285 6.80e-08 ***
## cargas      -0.168883   0.075703  -2.231   0.0300 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2646 on 52 degrees of freedom
## Multiple R-squared:  0.7444, Adjusted R-squared:  0.7248 
## F-statistic: 37.87 on 4 and 52 DF,  p-value: 8.003e-15

Sin embargo, cuando la variable dependiente es discreta (en este caso, binaria/dicotómica/dummy) no es adecuado estimar por MCO ya que se suelen presentar algunos problemas como los siguientes.

Falta de normalidad

Al representar los residuos se puede observar el comportamiento de una binomial, aunque en este caso no es claro:

residuos = mlp$residuals 
plot(residuos, xlab="Individuo", ylab="Residuos", col="red", lwd=2)
abline(h=0)

Representando de forma conjunta la variable dependiente y los residuos se observa que se distorsiona la distribución binomial incial:

plot(ts(cbind(credito, residuos)), plot.type="single", type="p", xlab="Individuo", ylab="Crédito (azul) y residuos (rojo)", col=c("blue","red"),lwd=2)

De igual forma, haciendo el contraste de normalidad de los residuos se suele rechazar la hipótesis nula de normalidad (p-valor<0.05), si bien, en este caso no se incumple esta hipótesis (p-valores por encima de 0.05):

shapiro.test(residuos) # H0: normalidad
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.96017, p-value = 0.05795
# install.packages("nortest") # instalarlo sólo una vez
library(nortest)
lillie.test(residuos) # H0: normalidad
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuos
## D = 0.1014, p-value = 0.1532

Falta de homocedasticidad

En la representación gráfica de los residuos en función de los ingresos se observa una tendencia creciente, lo cual nos hace pensar en la posible existencia de heterocedasticidad.

plot(ingresos,residuos^2,col="orange",lwd=3, xlab="Ingresos", ylab="Residuos al cuadrado")

plot(ingresos,residuos,col="brown",lwd=2, xlab="Ingresos", ylab="Residuos")

plot(laboral,residuos^2,col="orange",lwd=1, xlab="Situación laboral", ylab="Residuos al cuadrado") 

plot(laboral,residuos,col="brown",lwd=4, xlab="Situación laboral", ylab="Residuos") 

plot(cargas,residuos^2,col="orange",lwd=3, xlab="Cargas familiares", ylab="Residuos al cuadrado")

plot(cargas,residuos,col="brown",lwd=2, xlab="Cargas familiares", ylab="Residuos")

Esta sospecha se confirma al realizar el contraste de Goldfeld-Quandt (p-valor por debajo de 0.05, por lo que rechazo hipótesis nula) aunque no con el de Breusch-Pagan (p-valor superior a 0.05):

# install.packages("lmtest") # sólo hay que instalarlo una vez
library(lmtest) 
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(mlp) # H0: homocedasticidad
## 
##  studentized Breusch-Pagan test
## 
## data:  mlp
## BP = 8.2043, df = 4, p-value = 0.08437
gqtest(mlp, fraction=1/3, order.by=ingresos) # no me fio mucho
## 
##  Goldfeld-Quandt test
## 
## data:  mlp
## GQ = 2.0553e+31, df1 = 14, df2 = 14, p-value < 2.2e-16
## alternative hypothesis: variance increases from segment 1 to 2

Estimaciones fuera del intervalo [0, 1]

En el siguiente gráfico se observa que se obtienen estimaciones mayores que uno y menores que cero. Adviértase que en este caso se estima la probabilidad de devolver un crédito (en general, la probabilidad de que Y=1, siendo Y la variable dummy dependiente), por lo que estas estimaciones no tienen sentido:

estimaciones1=fitted.values(mlp) 
plot(estimaciones1, main="¿Datos estimados fuera de [0, 1]?", xlab="Individuo", ylab="Estimación de CREDITO", col="blue", lwd=2)
abline(a=1,b=0,col="red",lwd=2) # línea y=1
abline(a=0,b=0,col="red",lwd=2) # línea y=0

max(estimaciones1)
## [1] 1.273643
min(estimaciones1)
## [1] -0.1120744

Falta de representatividad del \(R^{2}\)

Observando la siguiente gráfica donde se representan los valores observados y estimados para la variable dependiente, se observa que la interpretación del coeficiente de determinación carece de sentido:

plot(ts(cbind(credito, estimaciones1)), plot.type="single", type="p", xlab="Individuo", ylab="Estimación de CREDITO (rojo) y CREDITO (azul)", col=c("blue","red"), lwd=2)

Para resolver estos problemas recurrimos a otro tipo de regresiones.

Regresión LOGIT

Un caso particular (que aparece con asiduidad en Economía) es aquel en el que la variable dependiente es binaria, es decir, toma únicamente dos valores. En tal caso, como hemos visto, la estimación por MCO no es adecuada, siendo necesaria una estimación logit o probit. Para ello recurrimos al código glm (esta función también se puede usar para corregir heteroscedasticidad y autocorrelación):

logit = glm(credito~ingresos+as.factor(laboral)+cargas, family=binomial("logit"))
summary(logit)
## 
## Call:
## glm(formula = credito ~ ingresos + as.factor(laboral) + cargas, 
##     family = binomial("logit"))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.51049  -0.00008   0.02909   0.25550   1.20168  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)  
## (Intercept)          -20.9772  4225.1521  -0.005   0.9960  
## ingresos               0.4847     0.2120   2.286   0.0223 *
## as.factor(laboral)1   18.2687  4225.1522   0.004   0.9966  
## as.factor(laboral)2   21.9597  4225.1523   0.005   0.9959  
## cargas                -2.1951     1.3133  -1.672   0.0946 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 79.001  on 56  degrees of freedom
## Residual deviance: 18.366  on 52  degrees of freedom
## AIC: 28.366
## 
## Number of Fisher Scoring iterations: 19

Ahora, por ejemplo, se tiene garantizado obtener estimaciones en el intervalo [0, 1]:

estimaciones2 = logit$fitted.values
plot(estimaciones2,main="¿Datos estimados fuera de [0, 1]?", xlab="Individuo", ylab="Estimación de CREDITO", col="blue",lwd=2)
abline(a=1,b=0,col="red",lwd=2) # línea y=1
abline(a=0,b=0,col="red",lwd=2) # línea y=0

max(estimaciones2)
## [1] 0.9995771
min(estimaciones2)
## [1] 1.702466e-10

Interpretación de las estimaciones obtenidas

En este caso, al ser el logit un modelo no lineal no es fácil interpretar las estimaciones de los coeficientes y los efectos marginales. Lo que si está claro es que:

  • Si el efecto marginal tiene signo positivo, entonces aumentos de la correspondiente variable implica aumento en la probabilidad de que Y=1.
  • Si el efecto marginal tiene signo negativo, entonces aumentos de la correspondiente variable implica disminución en la probabilidad de que Y=1.
  • El efecto marginal tiene el signo de la estimación del coeficiente.

En tal caso, puesto que la única variable con coeficiente significativamente distinto de cero es ingresos (p-valor=0.0223), y al ser este positivo, se tiene que conforme aumentan los ingresos mayor probabilidad de devolver el crédito.

Que el resto de coeficientes no sean significativamente distintos de cero indica que la probabilidad de dovelver el crédito es la misma para clientes en paro, con contrato temporal o indefinido y para clientes con carga y sin carga familiar.

Cálculo de probabilidad de éxito

Para el individuo medio se calcularía la probabilidad de éxito como sigue:

X = cbind(ingresos, laboral1, laboral2, cargas)
individuo.medio = colMeans(X) # calcula la media de una matriz por columnas
individuo.medio = c(1, individuo.medio) # añado 1 correspondiente a la cte
individuo.medio
##            ingresos  laboral1  laboral2    cargas 
## 1.0000000 6.1701754 0.3684211 0.3684211 0.4561404
coeficientes_logit = logit$coef # coeficientes de la regresión logit
coeficientes_logit
##         (Intercept)            ingresos as.factor(laboral)1 
##         -20.9771547           0.4846649          18.2687270 
## as.factor(laboral)2              cargas 
##          21.9597036          -2.1951490
z0 = sum(coeficientes_logit * individuo.medio) # ¡ojo a cómo se hace el producto!
z0
## [1] -4.166983
p0 = exp(z0)/(1+exp(z0))
p0 # probabilidad de éxito individuo medio
## [1] 0.0152624
1-p0 # probabilidad de fracaso individuo medio
## [1] 0.9847376

Para el individuo medio se tiene una probabilidad de devolver el crédito del 1.52%. Sin embargo, este individuo medio no tiene sentido ya que se ha considerado que en un 36.84% tiene contrato temporal, en un 36.84% contrato indefinido y tiene cargas familiares en un 45.61%.

Trabajemos con un individuo con ingresos medios que tenga contrato temporal y no tenga cargas familiares:

individuo1 = c(1, mean(ingresos), 1, 0, 0) 
individuo1
## [1] 1.000000 6.170175 1.000000 0.000000 0.000000
z1 = sum(coeficientes_logit * individuo1) 
z1
## [1] 0.2820397
p1 = exp(z1)/(1+exp(z1))
p1 # probabilidad de éxito individuo
## [1] 0.5700462
1-p1 # probabilidad de fracaso individuo 
## [1] 0.4299538

En este caso, la probabilidad de devolver el crédito para este individuo es del 57.004%.

Si se considerase que sí tiene cargas familiares la probabilidad pasaría a ser del 12.86%:

individuo2 = c(1, mean(ingresos), 1, 0, 1) 
individuo2
## [1] 1.000000 6.170175 1.000000 0.000000 1.000000
z2 = sum(coeficientes_logit * individuo2) 
z2
## [1] -1.913109
p2 = exp(z2)/(1+exp(z2))
p2 # probabilidad de éxito individuo 
## [1] 0.1286319
1-p2 # probabilidad de fracaso individuo 
## [1] 0.8713681

Cálculo e interpretación del odd

Si se divide la probabilidad de éxito entre la de fracaso se obtiene el número de veces más probable que ocurra el éxito que el fracaso (a esto se le conoce como odd). En este caso se estaría calculando el número de veces que es más probable que se devuelva el crédito frente a que no lo haga:

odd1 = p1/(1-p1)
odd1
## [1] 1.325831
odd2 = p2/(1-p2)
odd2
## [1] 0.1476207
1/odd2
## [1] 6.774119

Para el primer individuo, se obtiene que es 1.3258 veces más probable que devuelva el crédito frente a que no lo haga, mientras que para el segundo, es 6.774 veces menos probable que devuelva el crédito frente a que no lo haga.

En el caso del modelo logit, los odd también pueden calcularse como sigue:

exp(z1)
## [1] 1.325831
exp(z2)
## [1] 0.1476207
1/exp(z2)
## [1] 6.774119

Cálculo e interpretación del odd-ratio

Si se hace el cociente de dos odds, se obtiene el número de veces más probable que es devolver un crédito frente a no hacerlo asociado a los cambios producidos en las variables explicativas.

Así, con el cociente de los odds calculados anteriormente se obtiene que en un individuo con ingresos medios, con contrato temporal y con cargas (cargas=1) es 8.9813 veces menos probable que se devuelva el crédito frente a no hacerlo que cuando no tiene cargas (cargas=0):

odd_ratio=odd2/odd1
odd_ratio
## [1] 0.111342
1/odd_ratio
## [1] 8.981339

Adviertase que al obtener un odd-ratio inferior a 1 se ha cambiado la interpretación de más probable por menos probable.

¡Ojo! ¿Y si cambian los valores del resto de las variables? ¿Obtendremos el mismo odd-ratio? Puesto que el resto de variables permanecen constantes está claro que sí:

individuo3 = c(1, mean(ingresos), 0, 1, 0) 
individuo3
## [1] 1.000000 6.170175 0.000000 1.000000 0.000000
z3 = sum(coeficientes_logit * individuo3) 
z3
## [1] 3.973016
p3 = exp(z3)/(1+exp(z3))
p3 # probabilidad de éxito individuo
## [1] 0.9815309
1-p3 # probabilidad de fracaso individuo 
## [1] 0.01846907
odd3 = p3/(1-p3)
odd3
## [1] 53.14459
exp(z3)
## [1] 53.14459
individuo4 = c(1, mean(ingresos), 0, 1, 1) 
individuo4
## [1] 1.000000 6.170175 0.000000 1.000000 1.000000
z4 = sum(coeficientes_logit * individuo4) 
z4
## [1] 1.777867
p4 = exp(z4)/(1+exp(z4))
p4 # probabilidad de éxito individuo
## [1] 0.8554333
1-p4 # probabilidad de fracaso individuo 
## [1] 0.1445667
odd4 = p4/(1-p4)
odd4
## [1] 5.917223
exp(z4)
## [1] 5.917223
odd3/odd4
## [1] 8.981339

En el caso del modelo logit, el odd-ratio se puede obtener directamente calculando la exponencial de los coeficientes estimados:

coef = logit$coef[-1] # quito estimación de la cte
exp(coef) # odd-ratio última variable 
##            ingresos as.factor(laboral)1 as.factor(laboral)2 
##        1.623631e+00        8.590280e+07        3.443326e+09 
##              cargas 
##        1.113420e-01
1/exp(coef)[4]
##   cargas 
## 8.981339

Los odd-ratios anteriores se asocian a incrementos de una unidad en cada una variable. Si se quieren ver incrementos de más unidades hay que multiplicar el coeficiente estimado por dicha uniad y posteriormente calcular la exponencial:

exp(10*coef[1]) # odd-ratio ingresos, incremento 10 unidades
## ingresos 
##  127.313

En este caso, un individuo que tiene unos ingresos en 10000 euros superiores a otro (aumento de 10 unidades en la variable ingresos) es 127.313 veces más probable que devuelva el crédito frente a que no lo haga suponiendo que el resto de variables permanecen constantes.

Adviértase que al ser el resto de variables independientes binarias no tiene sentido considerar incrementos superiores a la unidad.

Efectos marginales

Aunque no sea fácil su interpretación cuantitativa (la cualitativa, el signo, ya se ha hecho al interpretar los coeficientes) vamos a calcular algunos efectos marginales.

Por ejemplo, para la variable cualitativa cargas habría que hacer las siguientes diferencias:

p2-p1
## [1] -0.4414143
p4-p3
## [1] -0.1260976

Puesto que hacen referencia a individuos distintos se obtienen valores distintos. En ambos casos se oberva que el pasar de no tener cargas familiares a sí hacerlo disminuye la probabilidad de devolver el crédito, cuestión que ya se sabía al haber obtenido un coeficiente estimado negativo para esta variable (este hecho no se ha interpretado anteriormente debido a que el coeficiente de la variable cargas no es significativamente distinto de cero).

Para la variable ingresos se tendría los siguienetes efectos marginales para cada uno de los individuos considerados hasta ahora:

em1 = p1*(1-p1)*coef[1]
em1
##  ingresos 
## 0.1187882
em2 = p2*(1-p2)*coef[1]
em2
##   ingresos 
## 0.05432403
em3 = p3*(1-p3)*coef[1]
em3
##    ingresos 
## 0.008785985
em4 = p4*(1-p4)*coef[1]
em4
##   ingresos 
## 0.05993713

Como se puede observar, dependen de los valores concretos de las variables explicativas usadas.

Bondad del ajuste

Como se ha comentado, el coeficiente de determinación deja de tener sentido. Es más, en la información anterior proporcionada por R no es calculado.

Una forma de estudiar la bondad del ajuste realizado puede ser la tabla de clasificación de aciertos y errores en la que se divide el número de casos clasificados correctamente entre el número de casos posibles.

Una función para calcularlos puede ser la siguiente:

tasa_aciertos<-function(corte,yest,y)
{
  A = 0
  B = 0
  C = 0
  D = 0
  n = length(y)
  for (i in 1:n){
    if ((yest[i]>=corte) & y[i]==1){A = A + 1}
    if ((yest[i]>=corte) & y[i]==0){B = B + 1}
    if ((yest[i]<corte) & y[i]==1){C = C + 1}
    if ((yest[i]<corte) & y[i]==0){D = D + 1}
  }
  tasa = ((A+D)/n)*100
  tasa
}
tasa_aciertos(mean(credito),fitted.values(logit),credito)
## [1] 92.98246

Para nuestro modelo se han clasificado correctamente el 92.98% de los casos.

Como corte se ha considerado la proporción de unos que hay en la variable dependiente (0.5087) en lugar de 0.5. Aunque en este caso prácticamente coinciden, de esta forma se adecúa el análisis a la información muestral en cada caso.

Finalmente, con el siguiente código se puede observar la evolución de la tasa de aciertos en función de diversos valores de corte:

j = 1
cortes = seq(0.01,0.99,0.01)
bondad = rep(0,length(cortes))
for (i in cortes)
     {
       bondad[j] = tasa_aciertos(i,fitted.values(logit),credito) 
       j = j + 1
}
plot(cortes,bondad,type="l",ylab="Tasa de casos clasificados correctamente", xlab="Corte", col="blue",lwd=2)

Para saber con detalle cómo se ha clasificado cada individuo:

estimacion = ifelse(estimaciones2 > mean(credito), 1, 0)
clasificacion = ifelse(credito == estimacion, "Correcto", "Incorrecto")
data.frame(credito, estimacion, clasificacion)
##    credito estimacion clasificacion
## 1        1          1      Correcto
## 2        1          1      Correcto
## 3        1          1      Correcto
## 4        1          1      Correcto
## 5        0          0      Correcto
## 6        0          0      Correcto
## 7        1          1      Correcto
## 8        1          1      Correcto
## 9        0          0      Correcto
## 10       0          0      Correcto
## 11       1          1      Correcto
## 12       0          0      Correcto
## 13       1          1      Correcto
## 14       1          1      Correcto
## 15       0          0      Correcto
## 16       0          0      Correcto
## 17       0          0      Correcto
## 18       0          0      Correcto
## 19       0          1    Incorrecto
## 20       1          1      Correcto
## 21       1          1      Correcto
## 22       1          0    Incorrecto
## 23       0          0      Correcto
## 24       0          0      Correcto
## 25       0          0      Correcto
## 26       1          1      Correcto
## 27       1          1      Correcto
## 28       1          1      Correcto
## 29       0          0      Correcto
## 30       1          1      Correcto
## 31       0          0      Correcto
## 32       1          1      Correcto
## 33       0          0      Correcto
## 34       1          1      Correcto
## 35       1          1      Correcto
## 36       0          0      Correcto
## 37       1          1      Correcto
## 38       1          0    Incorrecto
## 39       0          0      Correcto
## 40       0          0      Correcto
## 41       0          0      Correcto
## 42       1          1      Correcto
## 43       1          1      Correcto
## 44       1          1      Correcto
## 45       0          0      Correcto
## 46       0          0      Correcto
## 47       0          0      Correcto
## 48       1          1      Correcto
## 49       0          0      Correcto
## 50       1          1      Correcto
## 51       1          1      Correcto
## 52       0          0      Correcto
## 53       0          0      Correcto
## 54       0          0      Correcto
## 55       1          1      Correcto
## 56       1          1      Correcto
## 57       0          1    Incorrecto

Regresión PROBIT

Para la regresión probit hay que usar la siguiente orden:

probit = glm(credito~ingresos+as.factor(laboral)+cargas, family=binomial("probit"))
summary(probit)
## 
## Call:
## glm(formula = credito ~ ingresos + as.factor(laboral) + cargas, 
##     family = binomial("probit"))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.37719  -0.00010   0.00462   0.22918   1.28415  
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)  
## (Intercept)          -6.4991   570.0072  -0.011   0.9909  
## ingresos              0.2562     0.1016   2.521   0.0117 *
## as.factor(laboral)1   4.9853   570.0074   0.009   0.9930  
## as.factor(laboral)2   7.1625   570.0075   0.013   0.9900  
## cargas               -1.2033     0.6740  -1.785   0.0742 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 79.001  on 56  degrees of freedom
## Residual deviance: 18.330  on 52  degrees of freedom
## AIC: 28.33
## 
## Number of Fisher Scoring iterations: 18

Comparando las estimaciones de ambos modelos se observa que se obtienen el mismo signo en las estimaciones de los coeficientes, por lo que la interpretación del signo de los efectos marginales del probit coinciden con las del logit:

logit$coefficients
##         (Intercept)            ingresos as.factor(laboral)1 
##         -20.9771547           0.4846649          18.2687270 
## as.factor(laboral)2              cargas 
##          21.9597036          -2.1951490
probit$coefficients
##         (Intercept)            ingresos as.factor(laboral)1 
##           -6.499142            0.256224            4.985312 
## as.factor(laboral)2              cargas 
##            7.162486           -1.203317

Aunque la magnitud de las estimaciones de los coeficientes difieren, representando los valores estimados del modelo logit y probit se obtienen que prácticamente coinciden (se sitúan sobre la bisectriz del primer cuadrante):

plot(fitted.values(probit),fitted.values(logit),xlab="Valores ajustados modelo PROBIT",ylab="Valores ajustados modelo LOGIT",col="blue",lwd=2)
abline(a=0,b=1,col="red",lwd=2) # recta y=x

Interpretación coeficientes, odds, odds ratio y efectos marginales

La interpretación de estas estimaciones se hace como se ha indicado antes para el modelo logit, al igual que el cálculo e interpretación del odd y odd-ratio. Ahora bien, el uso de la exponencial para obtener estas medidas en este caso no es correcto ya que detrás del modelo probit está la función de distribución de la distribución normal y no la función logística (que depende exclusivamente de la exponencial).

Así, por ejemplo, para calcular la probabilidad de devolver un crédito para un individuo con ingresos medios que tenga contrato temporal y no tenga cargas familiares:

individuo1 = c(1, mean(ingresos), 1, 0, 0) 
individuo1
## [1] 1.000000 6.170175 1.000000 0.000000 0.000000
coeficientes_probit = probit$coefficients
coeficientes_probit
##         (Intercept)            ingresos as.factor(laboral)1 
##           -6.499142            0.256224            4.985312 
## as.factor(laboral)2              cargas 
##            7.162486           -1.203317
z1 = sum(coeficientes_probit * individuo1) 
z1
## [1] 0.06711711
p1 = pnorm(z1) # valor de la función distribición N(0,1)
p1 # probabilidad de éxito individuo 
## [1] 0.5267558
1-p1 # probabilidad de fracaso individuo 
## [1] 0.4732442
odd1 = p1/(1-p1) # odd
odd1
## [1] 1.113074

En este caso se tendría que la probabilidad de devolver el crédito del 52.67% y que es 1.113 veces más probable que se devuelva el crédito frente a que no lo haga. Mientras que si se considera que sí tiene cargas familiares, se tiene que la probabilidad de devolver el crédito frente a no hacerlo es del 12.79% y es 6.816 veces menos probable que se devuelva el crédito frente a que no se haga:

individuo2 = c(1, mean(ingresos), 1, 0, 1) 
individuo2
## [1] 1.000000 6.170175 1.000000 0.000000 1.000000
z2 = sum(coeficientes_probit * individuo2) 
z2
## [1] -1.1362
p2 = pnorm(z2) # valor de la función distribición N(0,1)
p2 # probabilidad de éxito individuo 
## [1] 0.1279365
1-p2 # probabilidad de fracaso individuo 
## [1] 0.8720635
odd2 = p2/(1-p2) # odd
odd2
## [1] 0.1467055
1/odd2
## [1] 6.816377

El odd-ratio vendría determinado por:

odd_ratio = odd2/odd1
odd_ratio
## [1] 0.1318021
1/odd_ratio
## [1] 7.58713

De forma que para este individuo (de ingresos medios y contrato temporal) si tiene cargas (cargas=1) familiares es 7.587 veces menos probable que devuelva el crédito frente a no hacerlo que cuando no las tiene (cargas=0).

Además, el efecto marginal se podría calcular como:

p2-p1
## [1] -0.3988192

Se tiene que en un individuo con ingresos medios y contrato temporal, la probabilidad de devolver el crédito disminuye un 39.88% cuando se pasa de no tener cargas familiares o tenerlas.

Finalmente, si se considera la variable cuantitativa refrerente a los ingresos, el efecto marginal se calcularía como:

em1 = dnorm(z1)*coeficientes_probit[2]  # dnorm() corresponde a la función de densidad de una N(0,1)
em1
##  ingresos 
## 0.1019886
em2 = dnorm(z2)*coeficientes_probit[2]
em2
##  ingresos 
## 0.0536048
em3 = dnorm(z3)*coeficientes_probit[2]
em3
##     ingresos 
## 3.818487e-05
em4 = dnorm(z4)*coeficientes_probit[2]
em4
##   ingresos 
## 0.02104593

Una vez más se observa que el efecto marginal depende de los valores usados en las variables explicativas.

Adviértase que se obtienen valores muy parecidos a los obtenidos con el modelo logit.

Bondad del ajuste

De igual forma, se puede obtener la tasa de aciertos para un umbral concreto (para la proporción de unos en la variable dependiente se obtiene una tasa de aciertos del 92.98%) y la evolución de la tasa de aciertos en función de diversos valores de corte:

tasa_aciertos(mean(credito),fitted.values(probit),credito)
## [1] 92.98246
j = 1
cortes = seq(0.01,0.99,0.01)
bondad = rep(0,length(cortes))
for (i in cortes)
     {
       bondad[j] = tasa_aciertos(i,fitted.values(probit),credito) 
       j = j + 1
}
plot(cortes,bondad,type="l",ylab="Tasa de casos clasificados correctamente", xlab="Corte", col="blue",lwd=2)

Una vez más se observa que los resultados son practicamente iguales a los del modelo logit.

Otra alternativa para medir la bondad del ajuste es la curva ROC, cuyo área bajo la curva puede interpretarse como la probabilidad de que ante un par de individuos, uno que devuelve el crédito y otro que no lo hace, la prueba los clasifique correctamente. Su valor máximo será el 1 y el mínimo 0.5. Un código para calcularla puede ser el siguiente:

CurvaROC<-function(corte,yest,y)
{
  A = 0
  B = 0
  C = 0
  D = 0
  n = length(y)
  for (i in 1:n){
    if ((yest[i]>=corte) & y[i]==1){A = A + 1}
    if ((yest[i]>=corte) & y[i]==0){B = B + 1}
    if ((yest[i]<corte) & y[i]==1){C = C + 1}
    if ((yest[i]<corte) & y[i]==0){D = D + 1}
  }
  especificidad = A/(A+B) # proporción de ceros clasificados correctamente
  sensibilidad = D/(C+D) # proporción de unos clasificados correctamente
  resultado = c(sensibilidad,1-especificidad)
  resultado
}

j = 1
cortes = seq(0.001,0.999,0.001)
sen = rep(0,length(cortes))
esp = rep(0,length(cortes))
for (i in cortes)
     {
       calculos = CurvaROC(i,fitted.values(probit),credito)
       sen[j] = calculos[1]
       esp[j] = calculos[2]
       j = j + 1
}
plot(esp,sen,type="l",ylab="Sensibilidad", xlab="1 - Especificidad", col="blue",lwd=2)

En este caso, simplemente hemos calculado la curva y no su área.

Finalmente, al igual que antes, para saber con detalle cómo se ha clasificado cada individuo:

estimacion = ifelse(fitted.values(probit) > mean(credito), 1, 0)
clasificacion = ifelse(credito == estimacion, "Correcto", "Incorrecto")
data.frame(credito, estimacion, clasificacion)
##    credito estimacion clasificacion
## 1        1          1      Correcto
## 2        1          1      Correcto
## 3        1          1      Correcto
## 4        1          1      Correcto
## 5        0          0      Correcto
## 6        0          0      Correcto
## 7        1          1      Correcto
## 8        1          1      Correcto
## 9        0          0      Correcto
## 10       0          0      Correcto
## 11       1          1      Correcto
## 12       0          0      Correcto
## 13       1          1      Correcto
## 14       1          1      Correcto
## 15       0          0      Correcto
## 16       0          0      Correcto
## 17       0          0      Correcto
## 18       0          0      Correcto
## 19       0          1    Incorrecto
## 20       1          1      Correcto
## 21       1          1      Correcto
## 22       1          0    Incorrecto
## 23       0          0      Correcto
## 24       0          0      Correcto
## 25       0          0      Correcto
## 26       1          1      Correcto
## 27       1          1      Correcto
## 28       1          1      Correcto
## 29       0          0      Correcto
## 30       1          1      Correcto
## 31       0          0      Correcto
## 32       1          1      Correcto
## 33       0          0      Correcto
## 34       1          1      Correcto
## 35       1          1      Correcto
## 36       0          0      Correcto
## 37       1          1      Correcto
## 38       1          0    Incorrecto
## 39       0          0      Correcto
## 40       0          0      Correcto
## 41       0          0      Correcto
## 42       1          1      Correcto
## 43       1          1      Correcto
## 44       1          1      Correcto
## 45       0          0      Correcto
## 46       0          0      Correcto
## 47       0          0      Correcto
## 48       1          1      Correcto
## 49       0          0      Correcto
## 50       1          1      Correcto
## 51       1          1      Correcto
## 52       0          0      Correcto
## 53       0          0      Correcto
## 54       0          0      Correcto
## 55       1          1      Correcto
## 56       1          1      Correcto
## 57       0          1    Incorrecto

Selección de modelos

Una vez descartada la estimación por MCO, ¿con cuál de los dos modelos quedarnos? ¿El logit o el probit?

Para responder a esta pregunta podríamos recurrir a los criterios de selección de modelos. En este caso, en los resultados ofrecidos por R se tiene el valor correspondiente al criterio de Akaike (AIC). Sin embargo, usando el comando mtable de la librería memisc se obtiene, entro otras cuestiones, el valor del criterio bayesiano de Schwarz (BIC):

# install.packages("memisc")
library(memisc)
## Loading required package: lattice
## Loading required package: MASS
## 
## Attaching package: 'memisc'
## The following objects are masked from 'package:stats':
## 
##     contr.sum, contr.treatment, contrasts
## The following object is masked from 'package:base':
## 
##     as.array
mtable(mlp, logit, probit)
## 
## Calls:
## mlp: lm(formula = credito ~ ingresos + as.factor(laboral) + cargas)
## logit: glm(formula = credito ~ ingresos + as.factor(laboral) + cargas, 
##     family = binomial("logit"))
## probit: glm(formula = credito ~ ingresos + as.factor(laboral) + cargas, 
##     family = binomial("probit"))
## 
## =============================================================
##                               mlp       logit      probit    
## -------------------------------------------------------------
##   (Intercept)              -0.006       -20.977     -6.499   
##                            (0.087)    (4225.152)  (570.007)  
##   ingresos                  0.045***      0.485*     0.256*  
##                            (0.009)       (0.212)    (0.102)  
##   as.factor(laboral): 1/0   0.203*       18.269      4.985   
##                            (0.099)    (4225.152)  (570.007)  
##   as.factor(laboral): 2/0   0.648***     21.960      7.162   
##                            (0.103)    (4225.152)  (570.008)  
##   cargas                   -0.169*       -2.195     -1.203   
##                            (0.076)       (1.313)    (0.674)  
## -------------------------------------------------------------
##   R-squared                 0.744                            
##   adj. R-squared            0.725                            
##   sigma                     0.265                            
##   F                        37.870                            
##   p                         0.000         0.000      0.000   
##   Log-likelihood           -2.478        -9.183     -9.165   
##   Deviance                  3.641        18.366     18.330   
##   AIC                      16.957        28.366     28.330   
##   BIC                      29.215        38.582     38.545   
##   N                        57            57         57       
##   Aldrich-Nelson R-sq.                    0.515      0.516   
##   McFadden R-sq.                          0.768      0.768   
##   Cox-Snell R-sq.                         0.655      0.655   
##   Nagelkerke R-sq.                        0.873      0.874   
##   phi                                     1.000      1.000   
##   Likelihood-ratio                       60.635     60.671   
## =============================================================

Para el logit se tiene valores iguales a 28.366 y 38.582, mientras que para el probit a 28.33 y 38.54. Por tanto, siguiendo estos criterios nos quedaríamos (por los pelos) con el probit.

Otra opción es elegir el modelo logit debido a la facilidad de cálculo que ofrece este modelo.

En cualquier caso, se ha observado que los resultados obtenidos por ambos modelos son muy similares, lo cual se ve reflejado en la obtención de un valor para el criterio AIC prácticamente idéntico.

Ejercicio propuesto

En la dirección web http://www.ugr.es/local/romansg/material/WebEco/04-Eco2/Ordenador/R/enfermedad.csv se tienen datos sobre las siguientes variables:

Usando las distintas herramientas presentadas, se pide: