Leyes anti-tabaco

Datos

Leemos la base de datos y la almacenamos en la variable smoke:

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

Además necsitaremos los siguientes paquetes de R:

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

Estadísticos descriptivos de las variables:

summary(smoke)
     educ          cigpric          white       
Min.   : 6.00   Min.   :44.00   Min.   :0.0000  
1st Qu.:10.00   1st Qu.:58.14   1st Qu.:1.0000  
Median :12.00   Median :61.05   Median :1.0000  
Mean   :12.47   Mean   :60.30   Mean   :0.8786  
3rd Qu.:13.50   3rd Qu.:63.18   3rd Qu.:1.0000  
Max.   :18.00   Max.   :70.13   Max.   :1.0000  
     age            income           cigs       
Min.   :17.00   Min.   :  500   Min.   : 0.000  
1st Qu.:28.00   1st Qu.:12500   1st Qu.: 0.000  
Median :38.00   Median :20000   Median : 0.000  
Mean   :41.24   Mean   :19305   Mean   : 8.686  
3rd Qu.:54.00   3rd Qu.:30000   3rd Qu.:20.000  
Max.   :88.00   Max.   :30000   Max.   :80.000  
   restaurn          smokes      
Min.   :0.0000   Min.   :0.0000  
1st Qu.:0.0000   1st Qu.:0.0000  
Median :0.0000   Median :0.0000  
Mean   :0.2466   Mean   :0.3841  
3rd Qu.:0.0000   3rd Qu.:1.0000  
Max.   :1.0000   Max.   :1.0000

Transformamos algunas variables para su uso posterior:

smoke <- within (smoke, {
  lincome <- log(income)
  lcigpric <- log(cigpric)
  age2 <- age^2
})

Modelo lineal de probabilidad

En primer lugar, estimamos un modelo lineal de probabilidad para smokes, una variable ficticia que indica si el individuo es fumador o no. Ahora sólo usaremos una variable explicativa: restaurn, una variable ficticia que indica si el individuo reside en un estado deonde está prohibido fumar en lugares públicos.

mlp <- lm (smokes ~ restaurn, data = smoke)
summary(mlp)
Call:
lm(formula = smokes ~ restaurn, data = smoke)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.4095 -0.4095 -0.3065  0.5905  0.6935 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.40954    0.01967  20.823  < 2e-16 ***
restaurn    -0.10301    0.03961  -2.601  0.00947 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.485 on 805 degrees of freedom
Multiple R-squared:  0.008332,	Adjusted R-squared:  0.007101 
F-statistic: 6.764 on 1 and 805 DF,  p-value: 0.009472

En un modelo lineal de probabilidad la varianza condicional del término de error cambia con las variables explicativas.

bptest(mlp, data = smoke)
    studentized Breusch-Pagan test

data:  mlp
BP = 9.1796, df = 1, p-value = 0.002447

Calculamos los errores típicos y los contrastes de hipótesis usando una matriz de covarianzas robusta a la heteroscedasticidad:

vrob <- vcovHC(mlp)
coeftest(mlp, vrob)
t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept)  0.409539   0.019976 20.5017 < 2.2e-16 ***
restaurn    -0.103007   0.038445 -2.6793  0.007528 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
waldtest(mlp, vcov = vrob)
Wald test

Model 1: smokes ~ restaurn
Model 2: smokes ~ 1
  Res.Df Df      F   Pr(>F)   
1    805                      
2    806 -1 7.1787 0.007528 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Otra alternativa es utilizar mínimos cuadrados ponderados. En un modelo lineal de probabilidad:

$$ E(u^2_i|x_i) = p_i (1 - p_i) $$

donde $$ p_i = \beta_0 + \beta_1 x_{1i} + \dots + \beta_k x_{ki} $$ Podemos estimar \(p_i\) con las predicciones del modelo lineal de probabilidad, reemplazando los parámetros por sus estimaciones en la expresión anterior, y utilizar como ponderaciones del método de MCP la inversa de

$$ h_i = \sqrt{\hat{p}_i (1 - \hat{p}_i)} $$

El procedimiento no funciona si la predicción para alguna observación es superior a 1 o inferior a 0.

phat <- fitted(mlp)
hhat <- sqrt(phat * (1 -phat))
mlp.mcp <- lm (smokes ~ restaurn, data = smoke, weights = 1/hhat)
summary(mlp.mcp)
Call:
lm(formula = smokes ~ restaurn, data = smoke, weights = 1/hhat)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-0.5840 -0.5840 -0.4514  0.8420  1.0213 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.40954    0.01981  20.670  < 2e-16 ***
restaurn    -0.10301    0.03895  -2.645  0.00834 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6967 on 805 degrees of freedom
Multiple R-squared:  0.008613,	Adjusted R-squared:  0.007381 
F-statistic: 6.994 on 1 and 805 DF,  p-value: 0.008339

Más variables explicativas

Las estimaciones anteriores indican un importante efecto de las leyes que restringen el uso del tabaco en espacios públicos. Ahora ampliamos el modelo para recoger otros factores que potencialmente influyen en el uso del tabaco:

  • lincome: logaritmo de la renta anual en dólares.
  • lcigpric: logaritmo del precio de un paquete de tabaco.
  • educ: años de educación.
  • white: variable ficticia que toma el valor 1 si el individuo no pertenece a una minoría racial.
  • age, age2: edad en años, edad al cuadrado.
mlp2 <- lm (smokes ~ restaurn + lincome + lcigpric + educ + white
            + age + age2 , data = smoke)
vrob <- vcovHC(mlp2)
coeftest(mlp2, vrob)
t test of coefficients:

               Estimate  Std. Error t value  Pr(>|t|)    
(Intercept)  0.65606243  0.87267620  0.7518 0.4524034    
restaurn    -0.10076163  0.03788692 -2.6595 0.0079816 ** 
lincome      0.01216204  0.02615171  0.4651 0.6420172    
lcigpric    -0.06892592  0.21002206 -0.3282 0.7428584    
educ        -0.02888022  0.00568094 -5.0837 4.613e-07 ***
white       -0.02568007  0.05129237 -0.5007 0.6167478    
age          0.01981580  0.00546916  3.6232 0.0003094 ***
age2        -0.00025979  0.00005784 -4.4915 8.116e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
waldtest(mlp2, vcov = vrob)
Wald test

Model 1: smokes ~ restaurn + lincome + lcigpric + educ + white + age + 
    age2
Model 2: smokes ~ 1
  Res.Df Df      F    Pr(>F)    
1    799                        
2    806 -7 10.385 1.627e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

En este caso no es posible utilizar MCP, dado que algunas predicciones caen fuera del intervalo (0, 1).

phat <- fitted (mlp2)
outside <- phat < 0 | phat > 1
sum(outside)
[1] 7
phat[outside]
          33          134          206          212          263 
-0.006829417 -0.019050044 -0.003965469 -0.264455016 -0.025876804 
         543          671 
-0.016104050 -0.035272189