Uso de una tarjeta de crédito

Datos

La base de datos ccard puede cargarse con la siguiente instrucción:

ccard <- read.csv2("http://jcpernias.com/ec1027/data/ccard.csv")

Además necsitaremos los siguientes paquetes de R:

ggplot2
gráficos.
sandwich
matrices de covarianzas robustas a heteroscedasticidad.
lmtest
contrastes de hipótesis.
library(ggplot2)
library(sandwich)
library(lmtest)

La base de datos ccard contiene información acerca del gasto mensual medio realizado con una tarjeta de crédito (variable expend, en dólares), la renta anual (variable income, en miles de dólares) y la edad (variable age, en años) de 72 clientes de un banco.

summary(ccard)
    expend             age            income      
Min.   :   9.58   Min.   :20.00   Min.   : 15.00  
1st Qu.:  67.60   1st Qu.:26.00   1st Qu.: 24.00  
Median : 158.32   Median :30.00   Median : 30.00  
Mean   : 262.53   Mean   :31.28   Mean   : 34.37  
3rd Qu.: 323.48   3rd Qu.:36.00   3rd Qu.: 39.70  
Max.   :1898.03   Max.   :55.00   Max.   :100.00

Parece existir una relación creciente entre los gastos y la renta. La dispersión también parece aumentar con la renta.

qplot(income, expend, data = ccard)

ccard-exp-inc.png

No está tan claro que exista una relación entre gastos y edad.

qplot(age, expend, data = ccard)

ccard-exp-age.png

Modelo de regresión

Planteamos el siguiente modelo de regresión donde la relación entre gastos y renta es cuadrática: $$expend_i = \beta_0 + \beta_1 age_i + \beta_2 income_i + \beta_3 income_i^2 + u_i$$

Creamos una nueva variable en la base de datos con los cuadrados de la renta anual:

ccard <- within (ccard, income2 <- income^2)

Estimamos los parámetros del modelo anterior por mínimos cuadrados ordinarios (MCO):

mco <- lm(expend ~ age + income + income2, data = ccard)
summary(mco)
Call:
lm(formula = expend ~ age + income + income2, data = ccard)

Residuals:
    Min      1Q  Median      3Q     Max 
-423.40 -126.27  -48.09   55.64 1444.35 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept) -255.40211  190.59372  -1.340  0.18469   
age           -2.51926    5.22156  -0.482  0.63102   
income        23.78550    7.91676   3.004  0.00372 **
income2       -0.15060    0.07418  -2.030  0.04626 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 282.9 on 68 degrees of freedom
Multiple R-squared:  0.2423,	Adjusted R-squared:  0.2089 
F-statistic: 7.248 on 3 and 68 DF,  p-value: 0.0002731

Heteroscedasticidad

Representamos gráficamente los residuos de la anterior estimación con respecto a la renta:

uhat <- resid(mco)
qplot(income, uhat, data=ccard)

ccard-uhat-inc.png

En el gráfico anterior la dispersión de los residuos es creciente con los niveles de renta. Para confirmar la presencia hetroscedasticidad llevamos a cabo el contraste de Breusch y Pagan con la función bptest del paquete lmtest:

bptest(mco)
    studentized Breusch-Pagan test

data:  mco
BP = 6.3361, df = 3, p-value = 0.09636

El p-value indica cierta evidencia en contra de la hipótesis de homoscedaticidad pero no podemos rechazar la hipótesis nula a un nivel de significación del 5%. Nuestro tamaño muestral es pequeño, \(n = 72\), y no está claro que la edad sea un determinate de los gastos. Utilicemos sólo \(income\) e \(income^2\) en la regresión auxiliar del contraste de Breusch y Pagan:

bptest(mco, ~ income + income2, data = ccard)
    studentized Breusch-Pagan test

data:  mco
BP = 6.2688, df = 2, p-value = 0.04353

En este caso, rechazamos la hipótesis de homoscedasticidad al \(5\%\).

Inferencia robusta a heteroscedasticidad

La presencia de heteroscedasticidad invalida las fórmulas usuales de los errores típicos de MCO y los contrastes \(F\). Para calcular una matriz de covarianzas válida aún en presencia de heteroscedasticidad usamos la función vcovHC del paquete sandwich:

vrob <- vcovHC(mco)
vrob
            (Intercept)           age        income      income2
(Intercept) 40011.55950 -383.24390717 -1473.8919874 12.016281056
age          -383.24391   11.79207186     1.6700109 -0.014378841
income      -1473.89199    1.67001092    74.7827288 -0.614609490
income2        12.01628   -0.01437884    -0.6146095  0.005110638

Los errores típicos de los estimadores son la raíz cuadrada de los elementos de la diagonal:

se.rob <- sqrt(diag(vrob))
se.rob
 (Intercept)          age       income      income2 
200.02889665   3.43395863   8.64770078   0.07148873

La comparación de los errores típicos robustos con los usuales no revela ningún patrón: en ocasiones los errores típicos robustos son mayores que los ordinarios y viceversa.

se.mco <- sqrt(diag(vcov(mco)))
cbind(mco = se.mco, rob = se.rob)
                     mco          rob
(Intercept) 190.59372201 200.02889665
age           5.22155907   3.43395863
income        7.91676408   8.64770078
income2       0.07418189   0.07148873

La función coeftest proporciona contrastes de significación individual similares a los obtenidos con summary, pero nos permite especificar un matriz de covarianzas robusta:

coeftest(mco, vrob)
t test of coefficients:

               Estimate  Std. Error t value Pr(>|t|)   
(Intercept) -255.402115  200.028897 -1.2768 0.206007   
age           -2.519265    3.433959 -0.7336 0.465694   
income        23.785501    8.647701  2.7505 0.007618 **
income2       -0.150595    0.071489 -2.1066 0.038848 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Para realizar el contraste de significación conjunta podemos usar la función waldtest del paquete lmtest:

waldtest(mco, vcov=vrob)
Wald test

Model 1: expend ~ age + income + income2
Model 2: expend ~ 1
  Res.Df Df      F    Pr(>F)    
1     68                        
2     71 -3 14.013 3.284e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Interpretación de los resultados

Nuestras estimaciones indican que la edad no tiene un efecto significativos sobre los gastos realizados con la tarjeta de crédito, una vez hemos tenido en cuenta la renta de los diferentes individuos. Por otro lado, existe una relación significativa entre los gastos y la renta. Pero esta relación es no lineal, lo que dificulta la interpretación de nuestras estimaciones.

En el siguiente gráfico se ha representado la relación entre la esperanza del gasto y la renta fijando la edad en su valor medio: $$ E(expend_i|income_i, \overline{age}) = \beta_0 + \beta_1 \overline{age} + \beta_2 income_i + \beta_3 income^2_i$$ La anterior expresión se estima reemplazando los parámetros \(\beta_0, \dots, \beta_3\) por sus estimaciones de MCO.

ccard-eexp-inc.png

El valor de la renta para el que se maximiza el gasto se puede obtener derivando la función de regresión con respecto a \(income\) e igualando a \(0\): $$ \frac{\partial E(expend)}{\partial income} = \beta_2 + 2 \beta_3 income^* = 0 $$ Por tanto: $$ income^* = -\frac{\beta_2}{2 \beta_3}$$ Evaluando esta última expresión con nuestras estimaciones:

b <- coef(mco)
income.star <- -b["income"] / (2 * b["income2"])
income.star
  income 
78.97159

Luego el efecto marginal de la renta sobre el gasto es negativo para rentas superiores a 78972 dólares. Sin embargo, este resultado ha de ser considerado con cierta cautela pues en nuestra muestra solo 2 de los 72 individuos tienen rentas anuales tan altas.

Otra forma de interpretar nuestros resultados es evaluar el efecto marginal de la renta para diferentes valores de \(income\):

income.values <- summary(ccard$income)
dexpend <- b["income"] + 2 * b["income2"] * income.values
cbind (income = income.values, "d expend/d income" = dexpend)
        income d expend/d income
Min.     15.00         19.267642
1st Qu.  24.00         16.556926
Median   30.00         14.749783
Mean     34.37         13.433580
3rd Qu.  39.70         11.828234
Max.    100.00         -6.333559

El efecto marginal de la renta sobre el gasto es decreciente con el nivel de renta. Para una renta anual de 24000 dólares, 1000 dólares adicionales de renta estań asociados con un incremento de los gastos mensuales de algo más de 16.5 dólares. El efecto marginal es de menos de 12 dólares mensuales para los clientes con una renta anual de 39700 dólares.

Mínimos cuadrados ponderados

Otra alternativa es usar un método de estimación eficiente en presencia de heteroscedasticidad. Para ello necesitamos un supuesto acerca de como se relaciona la varianza (condicional) del término de error con las variables explicativas. Supondremos que:

$$ E(u_i^2) = exp(\alpha_0 + \alpha_1 income_i + \alpha_2 income_i^2)$$

luhat2 <- log(uhat^2)
lh <- lm(luhat2 ~ income + income2, data = ccard)  
hhat <- exp(fitted(lh))

wls <- lm(expend ~ age + income + income2,
          weights = 1/hhat, data = ccard)
summary(wls)
Call:
lm(formula = expend ~ age + income + income2, data = ccard, weights = 1/hhat)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-2.5842 -1.2975 -0.4441  0.7374  8.2641 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)   
(Intercept) -137.97907   93.81174  -1.471  0.14596   
age           -0.55611    2.39763  -0.232  0.81728   
income        15.06579    4.46045   3.378  0.00121 **
income2       -0.08010    0.03664  -2.186  0.03226 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.08 on 68 degrees of freedom
Multiple R-squared:  0.8173,	Adjusted R-squared:  0.8093 
F-statistic: 101.4 on 3 and 68 DF,  p-value: < 2.2e-16

El siguiente gráfico compara los resultados de la estimación por MCP con la estimación por MCO. A pesar de las evidentes diferencias hay que tener en cuenta que en nuestra muestra hay muy pocas observaciones con una renta superior a 45000 dólares y en el rango de rentas que va desde 15000 a 45000 dólares las diferencias entre las dos estimaciones son menos acusadas.

ccard-eexp-inc-comp.png