IV sessions Week 9 Stat266 D Rogosa A. Observational data (adapted from Lab 3 Stat209) ################################################################################################ > mroz87 = read.table( "http://statweb.stanford.edu/~rag/stat209/Mroz87.dat", header = T) > names(mroz87) [1] "lfp" "hours" "kids5" "kids618" "age" "educ" [7] "wage" "repwage" "hushrs" "husage" "huseduc" "huswage" [13] "faminc" "mtr" "motheduc" "fatheduc" "unem" "city" [19] "exper" "nwifeinc" "wifecoll" "huscoll" > # Variable definition in lab script also linked, main page # lfp Dummy variable for labor-force participation. # hours Wife's hours of work in 1975. # kids5 Number of children 5 years old or younger. # kids618 Number of children 6 to 18 years old. # age Wife's age. # educ Wife's educational attainment, in years. # wage Wife's average hourly earnings, in 1975 dollars. # repwage Wife's wage reported at the time of the 1976 interview. # hushrs Husband's hours worked in 1975. # husage Husband's age. # huseduc Husband's educational attainment, in years. # huswage Husband's wage, in 1975 dollars. # faminc Family income, in 1975 dollars. # mtr Marginal tax rate facing the wife. # motheduc Wife's mother's educational attainment, in years. # fatheduc Wife's father's educational attainment, in years. # unem Unemployment rate in county of residence, in percentage points. # city Dummy variable = 1 if live in large city, else 0. # exper Actual years of wife's previous labor market experience. # nwifeinc Non-wife income. # wifecoll Dummy variable for wife's college attendance. # huscoll Dummy variable for husband's college attendance. > poswage = subset(mroz87, wage > 0) # my new data subset > poswage$logWage = log( poswage$wage ) # adding logWage to the data set for session > attach(poswage) > length(logWage) [1] 428 > table(lfp) #all the women in this subset are in the workforce > # onto fitting regression (predictor educ) > lm.posWage = lm(logWage ~ educ) > summary(lm.posWage) Call: lm(formula = logWage ~ educ) Residuals: Min 1Q Median 3Q Max -3.10256 -0.31473 0.06434 0.40081 2.10029 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -0.1852 0.1852 -1.000 0.318 educ 0.1086 0.0144 7.545 2.76e-13 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 0.68 on 426 degrees of freedom Multiple R-Squared: 0.1179, Adjusted R-squared: 0.1158 F-statistic: 56.93 on 1 and 426 DF, p-value: 2.761e-13 > # we get a highly significant slope (but not big Rsq), > # year increase in educ, fit increases 11 cents hourly wage > cor(logWage,educ)^2 # R-squared for the OLS equation [1] 0.1178826 > # now the IV fit using fatheduc as instrument (omitted vars concern) > cor.test(educ, fatheduc) Pearson's product-moment correlation data: educ and fatheduc t = 9.4255, df = 426, p-value < 2.2e-16 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: 0.3337579 0.4908623 sample estimates: cor 0.4154030 ##### Examples from Wooldridge, Introductory Econometrics ####### Chapters 15 and 16 ## stata results available from http://fmwww.bc.edu/gstat/examples/wooldridge/wooldridge15.html > install.packages("AER") > library(AER) > help(ivreg) > ivreg1 = ivreg(logWage ~ educ| fatheduc) > summary(ivreg1) Call: ivreg(formula = logWage ~ educ | fatheduc) Residuals: Min 1Q Median 3Q Max -3.0870 -0.3393 0.0525 0.4042 2.0677 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.44110 0.44610 0.989 0.323 educ 0.05917 0.03514 1.684 0.093 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.6894 on 426 degrees of freedom Multiple R-Squared: 0.09344, Adjusted R-squared: 0.09131 Wald test: 2.835 on 1 and 426 DF, p-value: 0.09294 > # matches tsls Task1 ####### use the diagnostics option of ivreg (recent) > ivreg1 = ivreg(logWage ~ educ| fatheduc) > summary(ivreg1, diagnostics = TRUE) Call: ivreg(formula = logWage ~ educ | fatheduc) Residuals: Min 1Q Median 3Q Max -3.0870 -0.3393 0.0525 0.4042 2.0677 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.44110 0.44610 0.989 0.3233 educ 0.05917 0.03514 1.684 0.0929 . Diagnostic tests: df1 df2 statistic p-value Weak instruments 1 426 88.84 <2e-16 *** Wu-Hausman 1 425 2.47 0.117 Sargan 0 NA NA NA --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.6894 on 426 degrees of freedom Multiple R-Squared: 0.09344, Adjusted R-squared: 0.09131 Wald test: 2.835 on 1 and 426 DF, p-value: 0.09294 > confint(ivreg1) # parameter confidence intervals as usual 2.5 % 97.5 % (Intercept) -0.433239996 1.3154468 educ -0.009703131 0.1280501 ##Investigate the 'Weak instruments' entry, matches cor.test > cor(educ, fatheduc) [1] 0.415403 ## just a test (like in the main lab of correlation between predictor and instrument # the Hausman test is often used to show there's a difference between OLS and IV, # but what does that tell you? # ivreg diagnostics finds discrepancy non-signicance even though IV est cuts the OLS estimate in half # and IV estimate non-significant ## to see the details of the Hausman test (and do it by hand) see pdf page 46 onward from http://personal.rhul.ac.uk/uhte/006/ec2203/Lecture%2015_IVestimation.pdf also the Basel class notes linked in week 6 ## The Sargan test is a statistical test used for testing over-identifying restrictions ## relevant to simultaneous eqs, which we investigated in identifiability in lab script # some package ivmodel display in week9 RQ ################################################################################################ ################################################################################################ B. Encouragement Design > library(foreign) > ses = read.dta("http://www.stat.columbia.edu/~gelman/arm/examples/sesame/sesame.dta") > names(ses) [1] "rownames" "id" "site" "sex" "age" "viewcat" "setting" "viewenc" "prebody" "prelet" "preform" "prenumb" "prerelat" "preclasf" "postbody" "postlet" "postform" [18] "postnumb" "postrelat" "postclasf" "peabody" "agecat" "encour" "_Isite_2" "_Isite_3" "_Isite_4" "_Isite_5" "regular" > dim(ses) [1] 240 28 > attach(ses) > table(encour) encour 0 1 88 152 > s1 = ivreg(postlet ~ as.numeric(viewcat) | encour, x = TRUE) > summary(s1, diagnostics = TRUE) Call: ivreg(formula = postlet ~ as.numeric(viewcat) | encour, x = TRUE) Residuals: Min 1Q Median 3Q Max -20.136 -8.830 -3.765 8.864 29.864 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 15.394 6.270 2.455 0.0148 * as.numeric(viewcat) 4.435 2.433 1.823 0.0695 . Diagnostic tests: df1 df2 statistic p-value Weak instruments 1 238 20.823 8.07e-06 *** Wu-Hausman 1 237 0.456 0.5 Sargan 0 NA NA NA --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 11.78 on 238 degrees of freedom Multiple R-Squared: 0.2282, Adjusted R-squared: 0.225 Wald test: 3.324 on 1 and 238 DF, p-value: 0.06954 > confint(s1) 2.5 % 97.5 % (Intercept) 3.1048997 27.683846 as.numeric(viewcat) -0.3328545 9.203703 > install.packages("ivpack") > library(ivpack) > robust.se(s1) [1] "Robust Standard Errors" t test of coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 15.3944 6.3210 2.4354 0.01561 * as.numeric(viewcat) 4.4354 2.4437 1.8150 0.07078 . --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 > anderson.rubin.ci(s1) $confidence.interval [1] "[ -1.28052883742684 , 9.37441745247359 ]" > ################################################################################################ ################################################################################################ C. Compliance Adjustments # binary compliance or measured compliance, week 9 RQ1 Artificial data in the image of Efron-Feldman # cholesterol reduction outcome measure > library(AER) > efdat = read.table(file = "http://web.stanford.edu/~rag/stat209/hw7efdata", header = T) > attach(efdat) > head(efdat) comp G Y 1 0.528 0 -9.57 2 0.862 1 58.50 3 0.980 0 10.10 4 0.673 1 33.90 5 0.660 0 3.60 6 0.551 1 29.10 > table(G) G 0 1 150 150 ###### got data > t.test(Y ~ G) # ITT estimate about 20, (not using compliance info) # approximating E-F paper analysis Welch Two Sample t-test data: Y by G t = -14.9477, df = 275.264, p-value < 2.2e-16 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -22.52955 -17.28584 sample estimates: mean in group 1 mean in group 2 9.663234 29.570932 # try an IV with compliance as a measured variable # and random assignment as instrument (analog to CACE) # Doesn't seem to work well, but see below > caceIV = ivreg(Y ~ comp|G) > confint(caceIV) 2.5 % 97.5 % (Intercept) -4262.740 2386.813 comp -3991.215 7220.521 # pretty wide > summary(caceIV) Call: ivreg(formula = Y ~ comp | G) Residuals: Min 1Q Median 3Q Max -634.86 -202.04 22.26 196.10 840.55 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -938 1696 -0.553 0.581 comp 1615 2860 0.565 0.573 Residual standard error: 305.4 on 298 degrees of freedom Multiple R-Squared: -399.6, Adjusted R-squared: -400.9 Wald test: 0.3187 on 1 and 298 DF, p-value: 0.5728 # but as you probably presumed that's not really the right analysis # the research question is, what is the effect of the drug? # the correct IV picture is as displayed in the Greenland handout # in the treatment group dose of the drug is confounded with compliance # in the control group dose of the drug is 0 (regardless of compliance level) # so the correct analog to AIR is > caceIV2 = ivreg(Y ~ I(comp*G)|G) > summary(caceIV2) Call: ivreg(formula = Y ~ I(comp * G) | G) Residuals: Min 1Q Median 3Q Max -28.5645 -6.1078 -0.2812 6.5362 29.6355 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 9.6645 0.8059 11.99 <2e-16 *** I(comp * G) 33.2189 1.9021 17.46 <2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 9.871 on 298 degrees of freedom Multiple R-Squared: 0.5814, Adjusted R-squared: 0.58 Wald test: 305 on 1 and 298 DF, p-value: < 2.2e-16 > # not far from the E-F paper estimate of 34.5 for ("CACE") effect (cholesterol reduction) for perfect compliers (these data are an imitation not exact recreation of E-F) > confint(caceIV2) 2.5 % 97.5 % (Intercept) 8.084849 11.24406 I(comp * G) 29.490842 36.94686 > # nice, tight interval ################################################ # Now go back to conventional Rubin approach with compliance as 0,1 > efdat$compT = efdat$comp > .8 # comp indicator # calculate simpler form of CACE with traditional .8 cutpoint, # yields with low compliance, 15%, maybe not the best choice of cut > mean(efdat$compT) [1] 0.15 > 20/.15 # divide ITT by proportion defined as compliant [1] 133.3333 > # wow > caceIVT = ivreg(Y ~ compT|G, data = efdat) > summary(caceIVT) Call: ivreg(formula = Y ~ compT | G, data = efdat) Residuals: Min 1Q Median 3Q Max -872.8 131.8 143.8 155.9 185.8 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -129.7 306.5 -0.423 0.673 compTTRUE 995.3 2039.0 0.488 0.626 # good bit diff than ITT of 20, CACE above is 133 Residual standard error: 353.2 on 298 degrees of freedom Multiple R-Squared: -534.8, Adjusted R-squared: -536.6 Wald test: 0.2383 on 1 and 298 DF, p-value: 0.6258 > confint(caceIVT) 2.5 % 97.5 % (Intercept) -730.4593 471.1126 compTTRUE -3001.0986 4991.6386 # here shows CACE can be anywhere for IV version # but as in above the right IV analysis uses dose (which equals compT in treat, 0 in control) > efdat$compT = efdat$comp > .8 # comp indicator > attach(efdat) > caceIV2 = ivreg(Y ~ I(compT*G)|G) # pretending dose is (0,1) > summary(caceIV2) # gives about usual case, wayoff from E-F effect for perfect compliance Call: ivreg(formula = Y ~ I(compT * G) | G) Residuals: Min 1Q Median 3Q Max -114.173 -4.502 5.436 16.561 46.436 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 9.664 2.447 3.950 9.78e-05 *** I(compT * G) 124.409 21.628 5.752 2.18e-08 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 29.97 on 298 degrees of freedom Multiple R-Squared: -2.858, Adjusted R-squared: -2.871 Wald test: 33.09 on 1 and 298 DF, p-value: 2.184e-08 > confint(caceIV2) 2.5 % 97.5 % (Intercept) 4.868521 14.46039 I(compT * G) 82.018297 166.79920 > > cor(efdat$G, efdat$compT) [1] 0.0280056 # G a pretty weak instrument here > cor.test(efdat$G, as.numeric(efdat$compT)) Pearson's product-moment correlation data: efdat$G and as.numeric(efdat$compT) t = 0.4836, df = 298, p-value = 0.629 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.08550641 0.14079991 sample estimates: cor 0.0280056 > cor.test(efdat$G, efdat$comp) Pearson's product-moment correlation data: efdat$G and efdat$comp t = 0.5565, df = 298, p-value = 0.5783 alternative hypothesis: true correlation is not equal to 0 95 percent confidence interval: -0.08131855 0.14493103 sample estimates: cor 0.03221898