Peso de los recién nacidos y consumo de tabaco

Puedes descargar aquí el script de R con el que se han obtenido los siguientes resultados.

Datos

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

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

Además necsitaremos los siguientes paquetes de R:

lmtest
contrastes de hipótesis.
library(lmtest)

Modelo básico

En primer lugar transformamos bwght. Esta variable contiene el peso de los recién nacidos en onzas. En nuestros modelos de regresión la utilizaremos en logaritmos.

bwght$lbwght <- log(bwght$bwght)

Nuestro primer modelo de regresión explica el peso de los recién nacidos sin tener en cuenta el consumo de tabaco. Como variables explicativas usamos faminc, la renta familiar anual en miles de dólares, parity, el orden de nacimiento del bebé, white, que toma el valor 1 para los bebés blancos, y male que toma el valor 1 para los varones.

mod0 <- lm(lbwght ~ faminc + parity + white + male, data = bwght)
summary(mod0)
Call:
lm(formula = lbwght ~ faminc + parity + white + male, data = bwght)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.58253 -0.09029  0.02110  0.12387  0.84328 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4.6583428  0.0168890 275.821  < 2e-16 ***
faminc      0.0008069  0.0002835   2.846  0.00449 ** 
parity      0.0143522  0.0056726   2.530  0.01151 *  
white       0.0518492  0.0129070   4.017 6.21e-05 ***
male        0.0271767  0.0101134   2.687  0.00729 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.188 on 1383 degrees of freedom
Multiple R-squared:  0.03079,	Adjusted R-squared:  0.02798 
F-statistic: 10.98 on 4 and 1383 DF,  p-value: 9.057e-09

Contrastes de especificación:

bptest(mod0)
    studentized Breusch-Pagan test

data:  mod0
BP = 5.8541, df = 4, p-value = 0.2103
resettest(mod0)
    RESET test

data:  mod0
RESET = 0.0068, df1 = 2, df2 = 1381, p-value = 0.9932

Ni el contraste de heteroscedasticidad de Breusch-Pagan ni el contraste RESET de Ramsey indican problemas en la especificación de nuestro modelo básico.

La variable parity, de acuerdo con las estimaciones anteriores, tiene un efecto positivo y significativo al 5%. Sin embargo, el valor que toma implica que el peso de los recién nacidos aumenta un 1.4% con el orden de nacimiento (el segundo hijo pesa en términos medios un 1.4% más que el primero, el tercero un 1.4% más que el segundo, etc). Sin embargo, una idea ampliamente extendida es que sólo existen diferencias de peso con respecto al primer hijo. Si en nuestra muestra el valor de parity tomara como mucho el valor 2, nuestras estimaciones reflejarían la diferencia de peso entre el primer y el segundo hijo y sería compatible con la idea de que es el primer hijo el que pesa menos. Podemos examinar los valores que toma la variable parity con la función table:

table(bwght$parity)
  1   2   3   4   5   6 
795 389 146  39  15   4

Casi un 15% de los recién nacidos de nuestra muestra tiene un valor de parity superior a 2. Como alternativa a nuestro modelo inicial plantearemos un modelo donde reemplazamos la variable parity por la variable first, una variable ficticia que toma el valor 1 en los casos en que el recién nacido sea el primer hijo y 0 en caso contrario:

bwght$first <- ifelse(bwght$parity == 1, 1, 0)

Las nuevas estimaciones son:

mod1 <- lm(lbwght ~ faminc + first + white + male, data = bwght)
summary(mod1)
Call:
lm(formula = lbwght ~ faminc + first + white + male, data = bwght)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.57943 -0.08856  0.02025  0.12455  0.84659 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.6982057  0.0141734 331.480  < 2e-16 ***
faminc       0.0007821  0.0002830   2.764  0.00579 ** 
first       -0.0282397  0.0102287  -2.761  0.00584 ** 
white        0.0523175  0.0129080   4.053 5.34e-05 ***
male         0.0273565  0.0101100   2.706  0.00690 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1879 on 1383 degrees of freedom
Multiple R-squared:  0.03164,	Adjusted R-squared:  0.02884 
F-statistic:  11.3 on 4 and 1383 DF,  p-value: 5.065e-09

El ajuste del modelo ha mejorado ligeramente y las estimaciones de los parámetros de las variables que hemos mantenido en el modelo apenas cambian. Las nuevas estimaciones implican que los primeros hijos pesan al nacer alrededor de un 2.8% menos que sus hermanos menores.

Es posible que al reemplazar parity por first hayamos desperdiciado mucha información. ¿Existen también diferencias de peso con respecto al segundo de los hijos? Creamos una variable ficticia que tome el valor 1 para los recién nacidos con parity igual a 2:

bwght$second <- ifelse(bwght$parity == 2, 1, 0)

Añadimos second al modelo de regresión anterior:

mod2 <- lm(lbwght ~ faminc + first + second + white + male,
           data = bwght)
summary(mod2)
Call:
lm(formula = lbwght ~ faminc + first + second + white + male, 
    data = bwght)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.57924 -0.08996  0.02047  0.12355  0.84645 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.7043491  0.0173665 270.887  < 2e-16 ***
faminc       0.0007919  0.0002835   2.793  0.00529 ** 
first       -0.0348070  0.0148207  -2.349  0.01899 *  
second      -0.0099729  0.0162835  -0.612  0.54034    
white        0.0524716  0.0129134   4.063 5.11e-05 ***
male         0.0273743  0.0101123   2.707  0.00687 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1879 on 1382 degrees of freedom
Multiple R-squared:  0.0319,	Adjusted R-squared:  0.0284 
F-statistic: 9.108 on 5 and 1382 DF,  p-value: 1.554e-08

La variable second no tiene un efecto significativo sobre el peso de los bebés.

Seleccionamos mod1 como nuestra especificación básica

summary(mod1)
Call:
lm(formula = lbwght ~ faminc + first + white + male, data = bwght)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.57943 -0.08856  0.02025  0.12455  0.84659 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.6982057  0.0141734 331.480  < 2e-16 ***
faminc       0.0007821  0.0002830   2.764  0.00579 ** 
first       -0.0282397  0.0102287  -2.761  0.00584 ** 
white        0.0523175  0.0129080   4.053 5.34e-05 ***
male         0.0273565  0.0101100   2.706  0.00690 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1879 on 1383 degrees of freedom
Multiple R-squared:  0.03164,	Adjusted R-squared:  0.02884 
F-statistic:  11.3 on 4 and 1383 DF,  p-value: 5.065e-09

Los contrastes de especificación no indican problemas:

bptest(mod1)
    studentized Breusch-Pagan test

data:  mod1
BP = 6.3732, df = 4, p-value = 0.173
resettest(mod1)
    RESET test

data:  mod1
RESET = 0.0109, df1 = 2, df2 = 1381, p-value = 0.9892

Efecto del consumo de tabaco

mod3 <- lm(lbwght ~ cigs + faminc + first + white + male,
           data = bwght)
summary(mod3)
Call:
lm(formula = lbwght ~ cigs + faminc + first + white + male, data = bwght)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.58629 -0.09019  0.02111  0.12020  0.84176 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.7147559  0.0144218 326.918  < 2e-16 ***
cigs        -0.0043277  0.0008528  -5.075 4.41e-07 ***
faminc       0.0005291  0.0002849   1.857  0.06349 .  
first       -0.0318701  0.0101636  -3.136  0.00175 ** 
white        0.0549876  0.0128048   4.294 1.87e-05 ***
male         0.0269963  0.0100210   2.694  0.00715 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1862 on 1382 degrees of freedom
Multiple R-squared:  0.04935,	Adjusted R-squared:  0.04591 
F-statistic: 14.35 on 5 and 1382 DF,  p-value: 1.014e-13
cor(bwght$cigs, bwght$faminc)
[1] -0.1730449
quantile(bwght$faminc)
 0%  25%  50%  75% 100% 
0.5 14.5 27.5 37.5 65.0
bwght$loinc <- ifelse(bwght$faminc < 14, 1, 0)
bwght$hiinc <- ifelse(bwght$faminc > 38, 1, 0)
mod4 <- lm(lbwght ~ cigs + loinc + hiinc + first + white + male,
           data = bwght)
summary(mod4)
Call:
lm(formula = lbwght ~ cigs + loinc + hiinc + first + white + 
    male, data = bwght)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.59756 -0.09019  0.02013  0.12013  0.83560 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.7376653  0.0147388 321.441  < 2e-16 ***
cigs        -0.0042176  0.0008514  -4.954 8.18e-07 ***
loinc       -0.0301339  0.0129852  -2.321  0.02045 *  
hiinc        0.0109317  0.0123565   0.885  0.37648    
first       -0.0325895  0.0101648  -3.206  0.00138 ** 
white        0.0505110  0.0128913   3.918 9.36e-05 ***
male         0.0279769  0.0100204   2.792  0.00531 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.186 on 1381 degrees of freedom
Multiple R-squared:  0.0523,	Adjusted R-squared:  0.04818 
F-statistic:  12.7 on 6 and 1381 DF,  p-value: 5.37e-14
mod4 <- lm(lbwght ~ cigs + loinc + first + white + male,
           data = bwght)
summary(mod4)
Call:
lm(formula = lbwght ~ cigs + loinc + first + white + male, data = bwght)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.59982 -0.08835  0.01966  0.12044  0.84266 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.7404284  0.0144030 329.129  < 2e-16 ***
cigs        -0.0042967  0.0008466  -5.075 4.40e-07 ***
loinc       -0.0331104  0.0125408  -2.640  0.00838 ** 
first       -0.0328603  0.0101594  -3.234  0.00125 ** 
white        0.0518942  0.0127952   4.056 5.28e-05 ***
male         0.0277469  0.0100162   2.770  0.00568 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.186 on 1382 degrees of freedom
Multiple R-squared:  0.05176,	Adjusted R-squared:  0.04833 
F-statistic: 15.09 on 5 and 1382 DF,  p-value: 1.881e-14
bptest(mod4)
    studentized Breusch-Pagan test

data:  mod4
BP = 7.5217, df = 5, p-value = 0.1846
resettest(mod4)
    RESET test

data:  mod4
RESET = 3.9959, df1 = 2, df2 = 1380, p-value = 0.0186
table(bwght$cigs)
   0    1    2    3    4    5    6    7    8    9   10   12   15 
1176    3    4    7    9   19    6    4    5    1   55    5   19 
  20   30   40   46   50 
  62    5    6    1    1
bwght <- within(bwght, {
  cigs5 <- ifelse(cigs > 0 & cigs <= 5, 1, 0)
  cigs10 <- ifelse(cigs > 5 & cigs <= 10, 1, 0)
  cigs20 <- ifelse(cigs > 10 & cigs <= 20, 1, 0)
  cigsabove20 <- ifelse(cigs > 20, 1, 0)
})
mod5 <- lm(lbwght ~ cigs5 + cigs10 + cigs20 + cigsabove20 + 
           loinc +  first + white + male, data = bwght)
summary(mod5)
Call:
lm(formula = lbwght ~ cigs5 + cigs10 + cigs20 + cigsabove20 + 
    loinc + first + white + male, data = bwght)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.59968 -0.08958  0.01996  0.11983  0.84003 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.74149    0.01444 328.309  < 2e-16 ***
cigs5       -0.03817    0.02935  -1.301  0.19359    
cigs10      -0.05922    0.02284  -2.593  0.00962 ** 
cigs20      -0.09439    0.02098  -4.500 7.36e-06 ***
cigsabove20 -0.10216    0.05198  -1.966  0.04955 *  
loinc       -0.03159    0.01258  -2.511  0.01217 *  
first       -0.03287    0.01016  -3.235  0.00125 ** 
white        0.05347    0.01282   4.171 3.22e-05 ***
male         0.02655    0.01005   2.643  0.00832 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.186 on 1379 degrees of freedom
Multiple R-squared:  0.05423,	Adjusted R-squared:  0.04874 
F-statistic: 9.883 on 8 and 1379 DF,  p-value: 1.915e-13
bwght$cigsabove10 <- ifelse(bwght$cigs > 10, 1, 0)
mod5 <- lm(lbwght ~ cigs5 + cigs10 + cigsabove10 + 
           loinc +  first + white + male, data = bwght)
summary(mod5)
Call:
lm(formula = lbwght ~ cigs5 + cigs10 + cigsabove10 + loinc + 
    first + white + male, data = bwght)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.59960 -0.08960  0.02025  0.12013  0.84001 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.74144    0.01443 328.516  < 2e-16 ***
cigs5       -0.03817    0.02934  -1.301  0.19345    
cigs10      -0.05922    0.02283  -2.594  0.00960 ** 
cigsabove10 -0.09542    0.01966  -4.853 1.35e-06 ***
loinc       -0.03159    0.01258  -2.511  0.01214 *  
first       -0.03287    0.01016  -3.237  0.00124 ** 
white        0.05354    0.01280   4.183 3.06e-05 ***
male         0.02653    0.01004   2.642  0.00833 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1859 on 1380 degrees of freedom
Multiple R-squared:  0.05421,	Adjusted R-squared:  0.04942 
F-statistic:  11.3 on 7 and 1380 DF,  p-value: 5.474e-14
bptest(mod5)
    studentized Breusch-Pagan test

data:  mod5
BP = 8.3247, df = 7, p-value = 0.3048
resettest(mod5)
    RESET test

data:  mod5
RESET = 2.2398, df1 = 2, df2 = 1378, p-value = 0.1069