library(tidyverse)
#install package
#install.packages("ivpack")
#load package
library(ivpack)
Data from the National Longitudinal Survey of Young Men (NLSYM) that was used by Card (1995), see the description. Nous considlwage
- traitement la variable educ
- variable instrumentale (IV) la variable near 4 year college
On veut d
data(card.data)
glimpse(card.data)
Observations: 3,010
Variables: 35
$ id <int> 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 19, 20, 22, 25, 26, 29, 30, 3...
$ nearc2 <int> 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0...
$ nearc4 <int> 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
$ educ <int> 7, 12, 12, 11, 12, 12, 18, 14, 12, 12, 9, 12, 11, 11, 16, 14, 12, 14, 10, 12, 18, 18,...
$ age <int> 29, 27, 34, 27, 34, 26, 33, 29, 28, 29, 28, 26, 24, 30, 31, 24, 34, 29, 26, 32, 32, 3...
$ fatheduc <int> NA, 8, 14, 11, 8, 9, 14, 14, 12, 12, 11, 11, 11, 11, NA, 15, 12, NA, 8, 8, 12, NA, 5,...
$ motheduc <int> NA, 8, 12, 12, 7, 12, 14, 14, 12, 12, 12, 6, 6, 6, 8, 12, 8, 12, 8, 8, 13, 8, 14, 12,...
$ weight <dbl> 158413, 380166, 367470, 380166, 367470, 380166, 367470, 496635, 367772, 480445, 38016...
$ momdad14 <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
$ sinmom14 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ step14 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ reg661 <int> 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1...
$ reg662 <int> 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0...
$ reg663 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ reg664 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ reg665 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ reg666 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ reg667 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ reg668 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ reg669 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ south66 <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ black <int> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ smsa <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1...
$ south <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0...
$ smsa66 <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
$ wage <int> 548, 481, 721, 250, 729, 500, 565, 608, 425, 515, 225, 400, 417, 217, 894, 300, 346, ...
$ enroll <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0...
$ KWW <int> 15, 35, 42, 25, 34, 38, 41, 46, 32, 34, 29, 34, 22, 27, 43, 36, 40, 35, 24, 35, 43, 5...
$ IQ <int> NA, 93, 103, 88, 108, 85, 119, 108, 96, 97, 84, 89, 93, 74, 116, NA, 93, 100, 91, 88,...
$ married <int> 1, 1, 1, 1, 1, 1, 1, 1, 4, 1, 1, 1, 1, 1, 1, 6, 1, 1, 6, 1, 1, 1, 1, 1, 1, 1, 6, 5, 1...
$ libcrd14 <int> 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1...
$ exper <int> 16, 9, 16, 10, 16, 8, 9, 9, 10, 11, 13, 8, 7, 13, 9, 4, 16, 9, 10, 14, 8, 10, 11, 10,...
$ lwage <dbl> 6.306275, 6.175867, 6.580639, 5.521461, 6.591674, 6.214608, 6.336826, 6.410175, 6.052...
$ expersq <int> 256, 81, 256, 100, 256, 64, 81, 81, 100, 121, 169, 64, 49, 169, 81, 16, 256, 81, 100,...
$ region <dbl> 661, 661, 661, 662, 662, 662, 662, 662, 662, 662, 669, 662, 662, 662, 662, 662, 662, ...
ggplot(card.data,aes(educ)) + geom_histogram()
Warning messages:
1: In native_encode(options$fig.path) :
some characters may not work under the current locale
2: In native_encode(options$fig.path) :
some characters may not work under the current locale
3: In native_encode(options$fig.path) :
some characters may not work under the current locale
4: In native_encode(options$fig.path) :
some characters may not work under the current locale
5: In native_encode(options$fig.path) :
some characters may not work under the current locale
6: In native_encode(options$fig.path) :
some characters may not work under the current locale
card.data %>% group_by(nearc4) %>% summarise(mean(educ))
Warning messages:
1: In native_encode(options$fig.path) :
some characters may not work under the current locale
2: In native_encode(options$fig.path) :
some characters may not work under the current locale
3: In native_encode(options$fig.path) :
some characters may not work under the current locale
card.data = card.data %>% mutate(educ12 = educ>12)
Warning messages:
1: In native_encode(options$fig.path) :
some characters may not work under the current locale
2: In native_encode(options$fig.path) :
some characters may not work under the current locale
3: In native_encode(options$fig.path) :
some characters may not work under the current locale
4: In native_encode(options$fig.path) :
some characters may not work under the current locale
5: In native_encode(options$fig.path) :
some characters may not work under the current locale
6: In native_encode(options$fig.path) :
some characters may not work under the current locale
7: In native_encode(options$fig.path) :
some characters may not work under the current locale
8: In native_encode(options$fig.path) :
some characters may not work under the current locale
card.data %>% group_by(nearc4) %>% summarise(moy = mean(educ12))
propcomp = card.data %>% group_by(nearc4) %>% summarise(moy = mean(educ12)) %>% pull(moy) %*% c(-1,1)
propcomp
[,1]
[1,] 0.1219293
card.data %>% group_by(nearc4) %>% summarise(moy = mean(lwage))
Warning messages:
1: In native_encode(options$fig.path) :
some characters may not work under the current locale
2: In native_encode(options$fig.path) :
some characters may not work under the current locale
3: In native_encode(options$fig.path) :
some characters may not work under the current locale
4: In native_encode(options$fig.path) :
some characters may not work under the current locale
5: In native_encode(options$fig.path) :
some characters may not work under the current locale
6: In native_encode(options$fig.path) :
some characters may not work under the current locale
7: In native_encode(options$fig.path) :
some characters may not work under the current locale
8: In native_encode(options$fig.path) :
some characters may not work under the current locale
itt = card.data %>% group_by(nearc4) %>% summarise(moy = mean(lwage)) %>% pull(moy) %*% c(-1,1)
itt
[,1]
[1,] 0.1559075
CACE
[,1]
[1,] 1.278672
# Stage 1
s1<-lm(educ12~nearc4,card.data)
predtx <-predict(s1, type = "response")
table(predtx)
predtx
0.422152560083588 0.54408183146614
957 2053
#Stage 2
lm(card.data$lwage~predtx)
Call:
lm(formula = card.data$lwage ~ predtx)
Coefficients:
(Intercept) predtx
5.616 1.279
ivpack
ivmodel1=ivreg(lwage ~ educ12, ~ nearc4, x=TRUE, data=card.data)
summary(ivmodel1)
Call:
ivreg(formula = lwage ~ educ12 | nearc4, data = card.data, x = TRUE)
Residuals:
Min 1Q Median 3Q Max
-2.28920 -0.53653 -0.03694 0.58481 2.04111
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.6157 0.1133 49.568 < 2e-16 ***
educ12TRUE 1.2787 0.2228 5.739 1.05e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.6941 on 3008 degrees of freedom
Multiple R-Squared: -1.445, Adjusted R-squared: -1.446
Wald test: 32.94 on 1 and 3008 DF, p-value: 1.047e-08
robust.se(ivmodel1)
[1] "Robust Standard Errors"
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.61570 0.11172 50.2637 < 2.2e-16 ***
educ12TRUE 1.27867 0.22036 5.8026 7.21e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
exper
, reg661
, reg662
, reg663
, reg664
, reg665
, reg666
, reg667
et reg668
ivmodel2=ivreg(lwage ~ educ12 + exper + reg661 + reg662 + reg663 + reg664 +
reg665+ reg666 + reg667 + reg668,
~ nearc4 + exper + reg661+ reg662 + reg663 + reg664 +
reg665 + reg666 +reg667 + reg668,
x=TRUE, data=card.data)
There were 11 warnings (use warnings() to see them)
summary(ivmodel2)
Call:
ivreg(formula = lwage ~ educ12 + exper + reg661 + reg662 + reg663 +
reg664 + reg665 + reg666 + reg667 + reg668 | nearc4 + exper +
reg661 + reg662 + reg663 + reg664 + reg665 + reg666 + reg667 +
reg668, data = card.data, x = TRUE)
Residuals:
Min 1Q Median 3Q Max
-2.55402 -0.39886 -0.01746 0.41183 1.69659
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.950242 0.366645 13.501 < 2e-16 ***
educ12TRUE 1.203660 0.307557 3.914 9.29e-05 ***
exper 0.081452 0.019639 4.148 3.45e-05 ***
reg661 0.009228 0.072198 0.128 0.898303
reg662 0.092530 0.055866 1.656 0.097770 .
reg663 0.109208 0.053864 2.027 0.042703 *
reg664 -0.023602 0.061442 -0.384 0.700905
reg665 -0.111628 0.065858 -1.695 0.090184 .
reg666 -0.135042 0.070416 -1.918 0.055234 .
reg667 -0.088088 0.068126 -1.293 0.196108
reg668 -0.256740 0.074531 -3.445 0.000579 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5709 on 2999 degrees of freedom
Multiple R-Squared: -0.6493, Adjusted R-squared: -0.6548
Wald test: 17.45 on 10 and 2999 DF, p-value: < 2.2e-16
robust.se(ivmodel2)
[1] "Robust Standard Errors"
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.950242 0.363037 13.6357 < 2.2e-16 ***
educ12TRUE 1.203660 0.304549 3.9523 7.920e-05 ***
exper 0.081452 0.019441 4.1898 2.873e-05 ***
reg661 0.009228 0.072869 0.1266 0.8992345
reg662 0.092530 0.056105 1.6492 0.0992065 .
reg663 0.109208 0.054991 1.9859 0.0471314 *
reg664 -0.023602 0.062978 -0.3748 0.7078594
reg665 -0.111628 0.065945 -1.6927 0.0906084 .
reg666 -0.135042 0.070452 -1.9168 0.0553593 .
reg667 -0.088088 0.068514 -1.2857 0.1986517
reg668 -0.256740 0.071859 -3.5729 0.0003587 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1