【Lecture Notes】Political Science 207

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

Political Science 207Winter Quarter 2014

• Home
  
• CV
  
• Contact
  

Final Exam

• Final Exam
  

Homeworks

• Homework #4 | Homework #4 Data
• Homework #3 | Homework #3 Data
• Homework #2 | Homework #2 Data
• Homework #1 | Homework #1 Data
  

Data and R Code

• R code for GIS | Data (San Francisco Voting Precincts) | R code for travel routing | Data (R Workspace for travel routing)
• R code for DiD, matching, RD | Data (Earned Income Tax Credit) | Data (Bridge Suicides) | Data (Incumbency Advantage)
• R code for multilevel models | Data (1992 U.S. Election Survey)
• R code for multiple imputation | Data (Central Coast Survey data) | Data (international development data)
• R code for selection models | Data (conflict data)
• R code for BTSCS data | Data (conflict data)
• R code for survival models | Data (coalition duration)
• R code for TSCS data | Data (Alvarez, Garrett, and Lange data)
• R code for instrumental variables | Data (Irish campaign data)
• R code for panel data | Data (British Election Study panel) | Data (Central bank data)
  

Readings

• Week 9: Difference-in-Difference, Matching, and Regression Discontinuity
  • Card, David, and Krueger, Alan B. 1994. "Minimum Wages and Employment: A Case Study of the Fast-Food Industry in New Jersey and Pennsylvania". American Economic Review 84: 772–793.
  • Ho, Daniel E., Kosuke Imai, Gary King, and Elizabeth A. Stuart. 2007. "Matching as Nonparametric Preprocessing for Reducing Model Dependence in Parametric Causal Inference." Political Analysis 15:199-236.
  • Gerber, Elisabeth R., and Daniel J. Hopkins. 2011. "When Mayors Matter: Estimating the Impact of Mayoral Partisanship on City Policy." American Journal of Political Science 55:326-339.
    
• Week 8: Multilevel Models
  • Jones, Bradford S. 2009. "Multilevel Modeling." In the The Oxford Handbook of Political Methodology, J. Box-Steffensmeier, H. Brady, and D. Collier (eds). New York: Oxford University Press.
  • Steenbergen, Marco R., and Bradford S. Jones. 2002. "Modeling Multilevel Data Structures." American Journal of Political Science 46:218-237.
    
• Week 7: Multiple Imputation for Missing Data
  • King, Gary, James Honaker, Anne Joseph, and Kenneth Scheve, 2001. "Analyzing Incomplete Political Science Data: An Alternative Algorithm for Multiple Imputation." American Political Science Review 95:49-70
  • Honaker, James, and Gary King. 2010. "What to do About Missing Values in Time Series Cross-Section Data." American Journal of Political Science 54:561-581.
    
• Week 6: Selection Models
  • Berinksy, Adam. 1999. "The Two Faces of Public Opinion." American Journal of Political Science 43:1209-1230.
  • Bushway, Shaun, Brian D. Johnson, and Lee Ann Slocum. 2007. "Is the Magic Still There? The Use of the Heckman Two-Step Correction for Selection Bias in Criminology." Journal of Quantitative Criminology 23:151-178.
  • Dubin, Jeffrey A., and Douglas Rivers, 1990. "Selection Bias in Linear Regression, Logit, and Probit Models." Sociologial Methods and Research 18:360-390.
    
• Week 5: Binary TSCS Models
  • Green, Donald P., Soo Yeon Kim, and David H. Yoon. 2001. "Dirty Pool" International Organization 55:441–468.
  • Beck, Nathaniel, and Jonathan N. Katz. 2001. "Throwing Out the Baby with the Bath Water: A Comment on Green, Kim, and Yoon." International Organization 55:487–495.
  • Beck, Nathaniel, Jonathan Katz, and Richard Tucker. 1998. "Taking Time Seriously: Time-Series-Cross-Section Analysis with a Binary Dependent Variable." American Journal of Political Science 42:1260-1288.
  • Carter, David B., and Curtis S. Signorino. 2010. "Back to the Future: Modeling Time Dependence in Binary Data." Political Analysis 18:271-292.
    
• Week 4: Survival Models
  • Box-Steffensmeier, Janet M. and Bradford S. Jones. 1997. "Time is of the Essence: Event History Models in Political Science."American Journal of Political Science 45:972–988.
  • Beck, Nathaniel, Jonathan Katz, and Richard Tucker. 1998. "Taking Time Seriously: Time-Series-Cross-Section Analysis with a Binary Dependent Variable." American Journal of Political Science 42:1260-1288.
    
• Week 3: Nonspherical Errors in TSCS Data
  • Beck, Nathaniel, and Jonathan N. Katz, 1995. "What to Do (and Not to Do) with Time-Series Cross-Section Data." American Political Science Review 89:634-647.
  • Beck, Nathaniel, and Jonathan N. Katz, 1996. "Nuisance vs. Substance: Specifying and Estimating Time-Series-Cross-Section Models." Political Analysis 6:1-36.
  • Alvarez, R. Michael, Geoffrey Garrett, and Peter Lange, 1991. "Government Partisanship, Labor Organization, and Macroeconomic Performance." American Political Science Review 85:539-556.
  • Beck, Nathaniel, Jonathan N. Katz, R. Michael Alvarez, Geoffrey Garrett, and Peter Lange, 1993. "Government Partisanship, Labor Organization, and Macroeconomic Performance: A Corrigendum." American Political Science Review 87:945-948.
    
• Week 2: Instrumental Variables
  • Benoit, Kenneth, and Michael Marsh. 2008. "The Campaign Value of Incumbency: A New Solution to the Puzzle of Less Effective Incumbent Spending." American Journal of Political Science , 52:874-890.
  • Gerber, Alan S., and Donald P. Green. 2000. "The Effects of Canvassing, Telephone Calls, and Direct Mail on Voter Turnout: A Field Experiment." American Political Science Review 94:653-663.
  • Jacobson, Gary. 1978. "The Effects of Campaign Spending in Congressional Elections." American Political Science Review72:469-491.
  • Alvarez, R. Michael, and Garrett Glasgow. 1999. "Two-Stage Estimation of Non-Recursive Choice Models." Political Analysis8:147-165.
    
• Week 1: Panel Data
  • Knack, Stephen, 1995. "Does 'Motor Voter' Work? Evidence from State-Level Data." Journal of Politics 57:796-811.
  • Stimson, James A., 1985. "Regression in Space and Time: A Statistical Essay." American Journal of Political Science 29:914-947.
  • Worrall, John L. 2010. "A User-friendly Introduction to Panel Data Modeling." Journal of Criminal Justice Education 21:182-196.
    
  
  

All course readings are available for download through the library website, either through JSTOR or our subscription to the electronic version of the journal. All campus URLs have access to these resources; if you want to download articles from home, see this.Course Materials

• Syllabus
  

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏1 回帖

关键词：matching travel

1. ## GIS in R ##
   
2. 
   
3. install.packages("maps")
   
4. library(maps)
   
5. install.packages("mapdata")
   
6. library(mapdata)
   
7. 
   
8. map("usa", col="blue")
   
9. 
   
10. map("state", col="blue")
    
11. 
    
12. map("county", col="blue")
    
13. 
    
14. install.packages("mapproj")
    
15. library(mapproj)
    
16. 
    
17. map("usa", col="blue", proj="sp_albers", par=c(30,40))  # conic projection
    
18. map("usa", col="blue", proj="sp_mercator")  # cylindrical projection
    
19. 
    
20. map("worldHires", "Canada", xlim=c(-141, -53), ylim=c(40, 85), col="gray", fill=TRUE)
    
21. map("worldHires", "Canada", xlim=c(-140, -110), ylim=c(48, 60), col="gray", fill=TRUE)
    
22. 
    
23. 
    
24. install.packages("rgdal")
    
25. library(rgdal)
    
26. 
    
27. sf.precincts <- readOGR(".", "sfprecincts")
    
28. 
    
29. redtags <- read.csv("redtags.csv")
    
30. 
    
31. 
    
32. install.packages("sp")  # should already be installed, but just in case
    
33. library(sp)
    
34. 
    
35. 
    
36. redtag.points <- SpatialPoints(cbind(redtags$lon, redtags$lat))
    
37. proj4string(redtag.points) <- proj4string(sf.precincts)
    
38. 
    
39. 
    
40. plot(sf.precincts)
    
41. plot(redtag.points, add=TRUE, col="red", pch=20)
    
42. 
    
43. sf.precincts@data$id <- seq(1,nrow(sf.precincts@data),1)
    
44. 
    
45. redtag.overlay <- over(redtag.points, sf.precincts)
    
46. redtag.overlay2 <- as.matrix(cbind(rownames(table(redtag.overlay$precinct)), table(redtag.overlay$precinct)))
    
47. colnames(redtag.overlay2) <- c("precinct", "redtags")
    
48. 
    
49. precinct.redtags <- merge(sf.precincts@data, redtag.overlay2, all.x=TRUE, sort=FALSE)
    
50. precinct.redtags$redtags <- as.numeric(as.character(precinct.redtags$redtags))
    
51. table(precinct.redtags$redtags)
    
52. 
    
53. sf.precincts@data <- precinct.redtags[order(precinct.redtags$id),]
    
54. sf.precincts@data$redtags[is.na(sf.precincts@data$redtags)] <- 0
    
55. 
    
56. 
    
57. install.packages("RColorBrewer")
    
58. library(RColorBrewer)
    
59. 
    
60. 
    
61. # plotting red tagged buildings
    
62. 
    
63. colors <- brewer.pal(5, "Reds")
    
64. brks <- c(0,1,2,5,10,30)
    
65. brknames <- c("0","1","2-4","5-9","10+")
    
66. 
    
67. plot(sf.precincts, col=colors[findInterval(sf.precincts$redtags, brks, all.inside=TRUE)])
    
68. plot(redtag.points, add=TRUE, col="blue", pch=1)
    
69. 
    
70. 
    
71. pdf("redtag_plot.pdf")
    
72. plot(sf.precincts, col=colors[findInterval(sf.precincts$redtags, brks, all.inside=TRUE)])
    
73. legend("topright", legend=brknames, fill=colors)
    
74. title("Red Tagged Buildings By Voting Precinct")
    
75. dev.off()
    
76. 
    
77. 
    
78. # plotting voter turnout
    
79. 
    
80. sf.precincts@data$turnout87[sf.precincts@data$turnout87==1] <- NA # set parks to missing
    
81. colors2 <- brewer.pal(5, "Greens")
    
82. brks2 <- c(0,quantile(sf.precincts@data$turnout87, na.rm=TRUE))
    
83. brknames2 <- c(paste(brks2[1], brks2[2], sep="-"), paste(brks2[2], brks2[3], sep="-"), paste(brks2[3], brks2[4], sep="-"), paste(brks2[4], brks2[5], sep="-"), paste(brks2[5], brks2[6], sep="-"))
    
84. 
    
85. colors2.map <- colors2[findInterval(sf.precincts$turnout87, brks2, all.inside=TRUE)]
    
86. colors2.map[is.na(colors2.map)] <- "#666666"
    
87. 
    
88. pdf("turnout_plot.pdf")
    
89. plot(sf.precincts, col=colors2.map)
    
90. legend("topright", legend=brknames2, fill=colors2)
    
91. title("Voter Turnout By Voting Precinct")
    
92. dev.off()
    
93. 
    
94. 
    
95. # calculation of distances
    
96. 
    
97. install.packages("geosphere")
    
98. library(geosphere)
    
99. 
    
100. 
     
101. # a simple spatial autocorrelation factor as in the bridges paper
     
102. 
     
103. precinct.dists <- as.matrix(dist(coordinates(sf.precincts)))
     
104. precinct.dists.inv <- 1/precinct.dists
     
105. diag(precinct.dists.inv) <- 0
     
106. 
     
107. install.packages("ape")
     
108. library(ape)
     
109. 
     
110. Moran.I(sf.precincts@data$redtags, precinct.dists.inv, na.rm=TRUE)
     
111. 
     
112. sf.precincts@data$sp.auto <- (precinct.dists.inv/1000) %*% sf.precincts@data$turnout89
     
113. 
     
114. reg.model1 <- lm(turnout89 ~ redtags + sp.auto, data = sf.precincts@data)
     
115. summary(reg.model1)

复制代码

1. ##########################
   
2. ## 12/8/2013
   
3. 
   
4. 
   
5. library(igraph)
   
6. library(ANN)
   
7. library(rgeos)
   
8. library(rgdal)
   
9. library(FNN)
   
10. library(geosphere)
    
11. 
    
12. sf.streets <- readOGR("C:\\R_group\\travel_routing", "tl_2012_06075_roads")
    
13. 
    
14. sf.streets@data$kph <- 1
    
15. sf.streets@data$kph[sf.streets@data$MTFCC=="S1100"] <- 100  # 65 mph
    
16. sf.streets@data$kph[sf.streets@data$MTFCC=="S1200"] <- 65   # 45 mph
    
17. sf.streets@data$kph[sf.streets@data$MTFCC=="S1400"] <- 40   # 25 mph
    
18. sf.streets@data$kph[sf.streets@data$MTFCC=="S1630"] <- 30
    
19. sf.streets@data$kph[sf.streets@data$MTFCC=="C3062"] <- 25
    
20. sf.streets@data$kph[sf.streets@data$MTFCC=="S1730"] <- 15
    
21. 
    
22. edge.count <- gIntersection(sf.streets,sf.streets)  ## printing shows 35542 edges
    
23. 
    
24. 
    
25. edges <- matrix(NA,100000,6)
    
26. 
    
27. row.start <- 0
    
28. 
    
29. # This loop takes about 80 minutes to run on my home machine
    
30. for (i in 1:length(sf.streets)) {
    
31.   temp.int <- gIntersection(sf.streets, sf.streets[i,])
    
32.   temp.edges <- do.call(rbind, lapply(temp.int@lines[[1]]@Lines, function(ls) {
    
33.     as.vector(t(ls@coords))  }))
    
34.   
    
35.   temp.lengths <- (distHaversine(temp.edges[,1:2], temp.edges[,3:4])/(sf.streets@data$kph[i]*1000)) * 60  # travel time in minutes
    
36.   
    
37.   edges[(row.start+1):(row.start+nrow(temp.edges)), 1:4] <- temp.edges
    
38.   edges[(row.start+1):(row.start+nrow(temp.edges)), 5] <- temp.lengths
    
39.   edges[(row.start+1):(row.start+nrow(temp.edges)), 6] <- i
    
40.   
    
41.   row.start <- row.start + nrow(temp.edges)
    
42.    
    
43. }
    
44. 
    
45. edges <- unique(na.omit(edges))
    
46. 
    
47. # start from here
    
48. 
    
49. library(igraph)
    
50. library(ANN)
    
51. library(rgeos)
    
52. library(rgdal)
    
53. library(FNN)
    
54. library(geosphere)
    
55. 
    
56. 
    
57. 
    
58. froms <- paste(edges[,1], edges[,2])
    
59. tos <- paste(edges[,3], edges[,4])
    
60. 
    
61. graph <- graph.edgelist(cbind(froms, tos), directed = FALSE)
    
62. E(graph)$weight <- edges[,5]
    
63. 
    
64. xy <- do.call(rbind, strsplit(V(graph)$name, " "))
    
65. 
    
66. V(graph)$x <- as.numeric(xy[, 1])
    
67. V(graph)$y <- as.numeric(xy[, 2])
    
68. xyg <- cbind(V(graph)$x, V(graph)$y)
    
69. 
    
70. ## pick route points here
    
71. ## Coit Tower -122.405833, 37.8025 -- 37.8025, -122.405833
    
72. ## Golden Gate Bridge (mid-span) -122.478611, 37.819722 -- 37.819722, -122.478611
    
73. ## Candlestick Park -122.386111, 37.713611 -- 37.713611,-122.386111
    
74. ## SF State -122.479722, 37.723333 -- 37.723333,-122.479722
    
75. 
    
76. from <- cbind(-122.478611, 37.819722)
    
77. to <- cbind(-122.405833, 37.8025)
    
78. 
    
79. 
    
80. ifrom <- get.knnx(xyg, from, 1)$nn.index[1, 1]
    
81. ito <- get.knnx(xyg, to, 1)$nn.index[1, 1]
    
82. v.path <- get.shortest.paths(graph, ifrom, ito, output = "vpath")[[1]]
    
83. e.path <- get.shortest.paths(graph, ifrom, ito, output = "epath")[[1]]
    
84. 
    
85. 
    
86. 
    
87. route <- xyg[v.path,]
    
88. totaldist <- 0
    
89. cardinal.directions <- NULL
    
90. for (i in 2:nrow(route)) {
    
91.   tempdist <- distHaversine(route[i-1,], route[i,])
    
92.   totaldist <- totaldist + tempdist
    
93.   direction <- bearing(route[i-1, ], route[i, ])
    
94.   cardinal.directions <- rbind(cardinal.directions, direction)
    
95. }
    
96. distance.km <- totaldist/1000 # distance in kilometers
    
97. distance.miles <- totaldist/1609.34  # distance in miles
    
98. travel.time <- sum(E(graph)$weight[v.path])  # travel time
    
99. 
    
100. street.names <- as.character(sf.streets@data[edges[e.path,6],2])
     
101. street.names[is.na(street.names)] <- "Unknown Road"
     
102. 
     
103. repeated <- matrix(FALSE,length(street.names),1)
     
104. for (j in 2:length(street.names)) {
     
105.   repeated[j] <- (street.names[j] == street.names[j-1])
     
106.   
     
107. }
     
108. 
     
109. compass.points <- cardinal.directions
     
110. compass.points[cardinal.directions>337.5 | cardinal.directions<22.5] <- "N"
     
111. compass.points[cardinal.directions>22.5 & cardinal.directions<67.5] <- "NE"
     
112. compass.points[cardinal.directions>67.5 & cardinal.directions<112.5] <- "E"
     
113. compass.points[cardinal.directions>112.5 & cardinal.directions<157.5] <- "SE"
     
114. compass.points[cardinal.directions>157.5 & cardinal.directions<202.5] <- "S"
     
115. compass.points[cardinal.directions>202.5 & cardinal.directions<247.5] <- "SW"
     
116. compass.points[cardinal.directions>247.5 & cardinal.directions<292.5] <- "W"
     
117. compass.points[cardinal.directions>292.5 & cardinal.directions<337.5] <- "NW"
     
118. 
     
119. all.turns <- cbind(street.names, compass.points, repeated)
     
120. turn.by.turn <- subset(all.turns, repeated==FALSE)[,1:2]
     
121. rownames(turn.by.turn) <- NULL
     
122. colnames(turn.by.turn) <- c("Street Name", "Direction")
     
123. 
     
124. route.results <- cbind(travel.time, distance.miles)
     
125. colnames(route.results) <- c("Travel Time (Minutes)", "Distance (Miles)")
     
126. print(route.results)
     
127. print(turn.by.turn)
     
128. 
     
129. plot(sf.streets, col="gray")
     
130. lines(V(graph)[v.path]$x, V(graph)[v.path]$y, lwd = 2, col = "red")

复制代码

1. ## PS 207 class 9
   
2. ## DiD, matching, RD
   
3. 
   
4. # earned income tax credit data
   
5. # EITC = earned income tax credit, a tax credit for lower income families with children.  
   
6. # Did more women work after it was expanded in 1994?
   
7. # urate = state unemployment rate
   
8. # finc = family income
   
9. # earn = earned income
   
10. # unearn = unearned income
    
11. 
    
12. library(foreign)
    
13. eitc <- read.dta("eitc.dta")
    
14. 
    
15. # Create two dummy variables to indicate before/after and treatment/control groups.
    
16. 
    
17. # the EITC went into effect in the year 1994
    
18. eitc$post93 <- as.numeric(eitc$year >= 1994)
    
19. 
    
20. # The EITC only affects women with at least one child, so the treatment group is all women with children.
    
21. eitc$anykids <- as.numeric(eitc$children >= 1)
    
22. 
    
23. # First calculate treatment effect through simple algebra
    
24. 
    
25. # Compute the four data points needed in the DiD calculation:
    
26. a <- mean(eitc$work[eitc$post93==0 & eitc$anykids==0], na.rm=T)
    
27. b <- mean(eitc$work[eitc$post93==0 & eitc$anykids==1], na.rm=T)
    
28. c <- mean(eitc$work[eitc$post93==1 & eitc$anykids==0], na.rm=T)
    
29. d <- mean(eitc$work[eitc$post93==1 & eitc$anykids==1], na.rm=T)
    
30. 
    
31. # Compute the effect of the EITC on the employment of women with children:
    
32. (d-c)-(b-a)
    
33. 
    
34. ## Now a DiD regression model
    
35. 
    
36. did.model1 <- lm(work ~ post93*anykids, data = eitc)
    
37. summary(did.model1)
    
38. 
    
39. did.model2 <- lm(work ~ post93*anykids + nonwhite + age + I(age^2) + ed + finc + I(finc-earn), data = eitc)
    
40. summary(did.model2)
    
41. 
    
42. ## placebo model
    
43. ## demonstrate we don't get a treatment effect when picking a different year for the "treatment"
    
44. 
    
45. ## subset to pre-treatment years
    
46. eitc.sub <- eitc[eitc$year <= 1993,]
    
47. 
    
48. # Create a new "after treatment" dummy variable
    
49. eitc.sub$post91 <- as.numeric(eitc.sub$year >= 1992)
    
50. 
    
51. # Run a placebo regression where placebo treatment = post91*anykids
    
52. did.model3 <- lm(work ~ post91*anykids, data = eitc.sub)
    
53. summary(did.model3)
    
54. 
    
55. #########################
    
56. ## Matching techniques ##
    
57. #########################
    
58. 
    
59. bridges <- read.dta("bridges.dta",  convert.factors=FALSE)
    
60. 
    
61. install.packages("MatchIt")
    
62. library(MatchIt)
    
63. 
    
64. # exact matching (fails) #
    
65. 
    
66. bridge.match <- matchit(totbexp ~ spcorr_all  + whitepct2564m + amindpct2564m + avgunemp + avgurban90 + totpop2564m, data=MatchData2, method="exact")
    
67. 
    
68. 
    
69. 
    
70. # propensity score matching #
    
71. 
    
72. bridge.match <- matchit(totbexp ~ spcorr_all  + whitepct2564m + amindpct2564m + avgunemp + avgurban90 + totpop2564m, data=MatchData2, method="nearest")
    
73. 
    
74. summary(bridge.match)
    
75. 
    
76. plot(bridge.match)
    
77. 
    
78. plot(bridge.match, type="jitter") # interactive!
    
79. 
    
80. plot(bridge.match, type="hist")
    
81. 
    
82. m.data <- match.data(bridge.match)
    
83. 
    
84. nrow(m.data)
    
85. 
    
86. 
    
87. # testing for balance with t-tests
    
88. t.test(avgunemp ~ totbexp, data=bridges)
    
89. t.test(avgunemp ~ totbexp, data=m.data)
    
90. 
    
91. 
    
92. # mahalanobis distance matching #
    
93. 
    
94. bridge.match2 <- matchit(totbexp ~ spcorr_all  + whitepct2564m + amindpct2564m + avgunemp + avgurban90 + totpop2564m, data=MatchData2, method="nearest", distance="mahalanobis")
    
95. 
    
96. summary(bridge.match2)
    
97. 
    
98. plot(bridge.match2)
    
99. 
    
100. # jitter and histogram plots are only for propensity scores
     
101. 
     
102. m2.data <- match.data(bridge.match2)
     
103. 
     
104. nrow(m2.data)
     
105. 
     
106. 
     
107. # propensity score matching, 2 matches per treated observation #
     
108. 
     
109. bridge.match3 <- matchit(totbexp ~ spcorr_all  + whitepct2564m + amindpct2564m + avgunemp + avgurban90 + totpop2564m, data=MatchData2, method="nearest", ratio=2)
     
110. 
     
111. summary(bridge.match3)
     
112. 
     
113. plot(bridge.match3)
     
114. 
     
115. plot(bridge.match3, type="jitter")
     
116. 
     
117. plot(bridge.match3, type="hist")
     
118. 
     
119. m3.data <- match.data(bridge.match3)
     
120. 
     
121. nrow(m3.data)
     
122. 
     
123. 
     
124. # propensity score matching,  discarding observations outside support #
     
125. 
     
126. bridge.match4 <- matchit(totbexp ~ spcorr_all  + whitepct2564m + amindpct2564m + avgunemp + avgurban90 + totpop2564m, data=MatchData2, method="nearest", discard="both")
     
127. 
     
128. summary(bridge.match4)
     
129. 
     
130. plot(bridge.match4)
     
131. 
     
132. plot(bridge.match4, type="jitter")
     
133. 
     
134. plot(bridge.match4, type="hist")
     
135. 
     
136. m4.data <- match.data(bridge.match4)
     
137. 
     
138. nrow(m4.data)
     
139. 
     
140. 
     
141. #install.packages("Zelig")
     
142. #library(Zelig)
     
143. # Zelig crashing, so using my own code
     
144. 
     
145. # Negative binomial model with unmatched data #
     
146. 
     
147. nbmod1 <- glm.nb(jumps2564_m  ~ spcorr_all  + whitepct2564m + amindpct2564m + avgunemp + avgurban90 + totbexp + offset(log(totpop2564m)),  data=bridges)
     
148. summary(nbmod1)
     
149. 
     
150. # Negative binomial model with propensity score matched data #
     
151. 
     
152. nbmod2 <- glm.nb(jumps2564_m  ~ spcorr_all  + whitepct2564m + amindpct2564m + avgunemp + avgurban90 + totbexp + offset(log(totpop2564m)),  data=m.data)
     
153. summary(nbmod2)
     
154. 
     
155. # Negative binomial model with mahalanobis distance matched data #
     
156. 
     
157. nbmod3 <- glm.nb(jumps2564_m  ~ spcorr_all  + whitepct2564m + amindpct2564m + avgunemp + avgurban90 + totbexp + offset(log(totpop2564m)),  data=m2.data)
     
158. summary(nbmod3)
     
159. 
     
160. # Negative binomial model with propensity score matching, 2 controls per treated #
     
161. 
     
162. nbmod4 <- glm.nb(jumps2564_m  ~ spcorr_all  + whitepct2564m + amindpct2564m + avgunemp + avgurban90 + totbexp + offset(log(totpop2564m)),  data=m3.data)
     
163. summary(nbmod4)
     
164. 
     
165. # Negative binomial model with propensity score matching, discarding observations w/o common support #
     
166. 
     
167. nbmod5 <- glm.nb(jumps2564_m  ~ spcorr_all  + whitepct2564m + amindpct2564m + avgunemp + avgurban90 + totbexp + offset(log(totpop2564m)),  data=m4.data)
     
168. summary(nbmod5)
     
169. 
     
170. x.t0 <- cbind(t(colMeans(cbind(1, bridges$spcorr_all, bridges$whitepct2564m, bridges$amindpct2564m, bridges$avgunemp, bridges$avgurban90))), 0)
     
171. x.t1 <- cbind(t(colMeans(cbind(1, bridges$spcorr_all, bridges$whitepct2564m, bridges$amindpct2564m, bridges$avgunemp, bridges$avgurban90))), 1)
     
172. offset.var <- mean(log(bridges$totpop2564m))
     
173. 
     
174. mod1.treat0 <- exp(nbmod1$coef%*%t(x.t0) + offset.var)
     
175. mod1.treat1 <- exp(nbmod1$coef%*%t(x.t1) + offset.var)
     
176. att1 <- mod1.treat1 - mod1.treat0
     
177. mod1.treat1
     
178. mod1.treat0
     
179. att1
     
180. 
     
181. mod5.treat0 <- exp(nbmod5$coef%*%t(x.t0) + offset.var)
     
182. mod5.treat1 <- exp(nbmod5$coef%*%t(x.t1) + offset.var)
     
183. att5 <- mod5.treat1 - mod5.treat0
     
184. mod5.treat1
     
185. mod5.treat0
     
186. att5
     
187. 
     
188. 
     
189. 
     
190. 
     
191. ##############################
     
192. ## Regression Discontinuity ##
     
193. ##############################
     
194. 
     
195. # incumbency advantage
     
196. votedata <- read.dta("LMBdata.dta")
     
197. 
     
198. votedata$treatment <- (votedata$lagdemvoteshare >= 0.5)
     
199. 
     
200. rd.model1 <- lm(demvoteshare ~ lagdemvoteshare + treatment, data=votedata)
     
201. summary(rd.model1)
     
202. 
     
203. plot(demvoteshare ~ lagdemvoteshare, data=votedata)
     
204. curve(rd.model1$coefficients[1] + rd.model1$coefficients[2]*x  + rd.model1$coefficients[3]*(x>0.5), from=0, to = 0.5, col="red", lw=2, add=TRUE)
     
205. curve(rd.model1$coefficients[1] + rd.model1$coefficients[2]*x  + rd.model1$coefficients[3]*(x>0.5), from=0.501, to=1, col="red", lw=2, add=TRUE)
     
206. 
     
207. 
     
208. rd.model2 <- lm(demvoteshare ~ lagdemvoteshare*treatment, data=votedata)
     
209. summary(rd.model2)
     
210. 
     
211. plot(demvoteshare ~ lagdemvoteshare, data=votedata)
     
212. curve(rd.model2$coefficients[1] + rd.model2$coefficients[2]*x  + rd.model2$coefficients[3]*(x>0.5) + rd.model2$coefficients[4]*(x>0.5)*x, from=0, to = 0.5, col="red", lw=2, add=TRUE)
     
213. curve(rd.model2$coefficients[1] + rd.model2$coefficients[2]*x  + rd.model2$coefficients[3]*(x>0.5) + rd.model2$coefficients[4]*(x>0.5)*x, from=0.501, to = 1, col="red", lw=2, add=TRUE)
     
214. 
     
215. 
     
216. rd.model3 <- lm(demvoteshare ~ lagdemvoteshare*treatment + I(lagdemvoteshare^2)*treatment, data=votedata)
     
217. summary(rd.model3)
     
218. 
     
219. plot(demvoteshare ~ lagdemvoteshare, data=votedata)
     
220. curve(rd.model3$coefficients[1] + rd.model3$coefficients[2]*x  + rd.model3$coefficients[3]*(x>0.5) +  rd.model3$coefficients[4]*x*x + rd.model3$coefficients[5]*(x>0.5)*x + rd.model3$coefficients[6]*(x>0.5)*x*x, from=0, to = 0.5, col="red", lw=2, add=TRUE)
     
221. curve(rd.model3$coefficients[1] + rd.model3$coefficients[2]*x  + rd.model3$coefficients[3]*(x>0.5) +  rd.model3$coefficients[4]*x*x + rd.model3$coefficients[5]*(x>0.5)*x + rd.model3$coefficients[6]*(x>0.5)*x*x, from=0.501, to = 1, col="red", lw=2, add=TRUE)
     
222. 
     
223. 
     
224. install.packages("rdd")
     
225. library(rdd)
     
226. 
     
227. 
     
228. rd.model4 <- RDestimate(demvoteshare ~ lagdemvoteshare, cutpoint = 0.5, data=votedata)
     
229. summary(rd.model4)
     
230. plot(rd.model4)
     
231. 
     
232. # placebo test around median of lagged vote share
     
233. 
     
234. rd.model5 <- RDestimate(demvoteshare ~ lagdemvoteshare, cutpoint = median(votedata$lagdemvoteshare, na.rm=T), data=votedata)
     
235. summary(rd.model5)
     
236. 
     
237. # trying different bandwidths
     
238. 
     
239. rd.model4$bw
     
240. 
     
241. rd.model.05 <- RDestimate(demvoteshare ~ lagdemvoteshare, cutpoint = 0.5, data=votedata, bw=0.05)
     
242. summary(rd.model.05)
     
243. plot(rd.model.05)
     
244. 
     
245. rd.model.025 <- RDestimate(demvoteshare ~ lagdemvoteshare, cutpoint = 0.5, data=votedata, bw=0.025)
     
246. summary(rd.model.025)
     
247. plot(rd.model.025)
     
248. 
     
249. rd.model.01 <- RDestimate(demvoteshare ~ lagdemvoteshare, cutpoint = 0.5, data=votedata, bw=0.01)
     
250. summary(rd.model.01)
     
251. plot(rd.model.01)

复制代码

1. ## multilevel model code for R ##
   
2. 
   
3. install.packages("lme4")
   
4. library(lme4)
   
5. 
   
6. 
   
7. 
   
8. #### First example -- exam scores
   
9. # normexam = test scores
   
10. # school = school id
    
11. # standLRT = individual score on a different test
    
12. # schavg = average intake score
    
13. 
    
14. install.packages("mlmRev")
    
15. library(mlmRev)
    
16. 
    
17. data(Exam)
    
18. head(Exam)
    
19. 
    
20. lmer(normexam ~ 1 + (1 | school), data=Exam) # random intercept only
    
21. 
    
22. lmer(normexam ~ standLRT + (1 | school), data=Exam) # ranomd intercept plus fixed effect
    
23. 
    
24. lmer(normexam ~ standLRT + (standLRT | school), data=Exam) # random intercept, random slope
    
25. 
    
26. lmer(normexam ~ standLRT + schavg + (1 + standLRT | school), data=Exam) # random intercept, individual and group predictor
    
27. 
    
28. lmer(normexam ~ standLRT * schavg + (1 + standLRT | school), data=Exam) # random intercept, cross-level interaction
    
29. 
    
30. 
    
31. 
    
32. 
    
33. 
    
34. 
    
35. 
    
36. data(InstEval)
    
37. # s = student number
    
38. # d = professor number
    
39. # studage = student age (# semesters enrolled)
    
40. # lectage = lecture age (# semesters lecture was in the past)
    
41. # service = dummy variable for course taught outside department
    
42. # dept = department number
    
43. # y = student rating, higher is better
    
44. 
    
45. InstEval$studage <- as.numeric(as.character(InstEval$studage)) # convert from factor to number
    
46. head(InstEval)
    
47. 
    
48. ## linear regression
    
49. 
    
50. linmod1 <- lm(y ~ studage + service, data=InstEval)
    
51. summary(linmod1)
    
52. 
    
53. linmod2 <- lm(y ~ studage + service + as.factor(dept), data=InstEval)
    
54. summary(linmod2)
    
55. 
    
56. 
    
57. ## linear mixed model, random intercept
    
58. 
    
59. mixedmod1 <- lmer(y ~ studage + service + (1 | dept), data=InstEval)
    
60. summary(mixedmod1)
    
61. 
    
62. ## linear regression, interaction model
    
63. 
    
64. linmod3 <- lm(y ~ studage + service + (studage * dept), data=InstEval)
    
65. summary(linmod3)
    
66. 
    
67. ## linear mixed model, random intercept and random coefficient on student age by department
    
68. ## note studage enters twice -- otherwise mean effect assumed to be zero
    
69. 
    
70. mixedmod2 <- lmer(y ~ service + studage + (1 + studage | dept), data=InstEval)
    
71. summary(mixedmod2)
    
72. 
    
73. # examine random and fixed effects
    
74. fixef(mixedmod2)
    
75. ranef(mixedmod2)
    
76. 
    
77. # plot results
    
78. 
    
79. dotplot(ranef(mixedmod2))
    
80. 
    
81. randeffs <- unlist(ranef(mixedmod2))
    
82. reglines <- cbind((randeffs[1:14] + fixef(mixedmod2)[1]), (randeffs[15:28] + fixef(mixedmod2)[2]))
    
83. 
    
84. plot(y ~ studage, ylim = c(2,3.5), data=InstEval)
    
85. for (i in 1:nrow(reglines)) {
    
86.   curve(reglines[i,1] + reglines[i,2] * x, col=i, add=TRUE)
    
87. }
    
88. 
    
89. 
    
90. # compare multilevel to separate regressions -- shrinkage estimator
    
91. 
    
92. 
    
93. mixedmod3 <- lmer(y ~ studage + (1 + studage | dept), data=InstEval)
    
94. summary(mixedmod3)
    
95. randeffs3 <- unlist(ranef(mixedmod3))
    
96. reglines3 <- cbind((randeffs3[1:14] + fixef(mixedmod3)[1]), (randeffs3[15:28] + fixef(mixedmod3)[2]))
    
97. 
    
98. alpha.hat <- NULL
    
99. beta.hat <- NULL
    
100. for(i in levels(InstEval$dept)){
     
101.   unit.lm <- lm(y ~ studage, data = subset(InstEval, dept == i) )
     
102.   alpha.hat <- append(alpha.hat, coef(unit.lm)[1])
     
103.   beta.hat <- append(beta.hat, coef(unit.lm)[2])
     
104. }
     
105. 
     
106. pooled.model <- lm(y ~ studage, data=InstEval)
     
107. pooled.beta <- coef(pooled.model)[2]
     
108. 
     
109. plot(reglines3[,2] ~ beta.hat)
     
110. curve(1*x, add=TRUE)
     
111. abline(v=pooled.beta)
     
112. 
     
113. 
     
114. 
     
115. ## Using Zelig
     
116. 
     
117. install.packages("ZeligMultilevel")
     
118. library(ZeligMultilevel)
     
119. 
     
120. 
     
121. # example 1
     
122. 
     
123. z.reg1 <- zelig(y ~ studage + service + tag(1 | dept), data=InstEval, model="ls.mixed")
     
124. summary(z.reg1)
     
125. 
     
126. x.high <- setx(z.reg1, studage=8)
     
127. x.low <- setx(z.reg1, studage=2)
     
128. 
     
129. s.reg1 <- sim(z.reg1, x=x.high, x1=x.low)
     
130. summary(s.reg1)
     
131. 
     
132. # example 2 -- note assumed intercept in random component
     
133. 
     
134. z.reg2 <- zelig(y ~ studage + service + tag(studage | dept), data=InstEval, model="ls.mixed")
     
135. summary(z.reg2)
     
136. 
     
137. x.high <- setx(z.reg2, studage=8)
     
138. x.low <- setx(z.reg2, studage=2)
     
139. 
     
140. s.reg2 <- sim(z.reg2, x=x.high, x1=x.low)
     
141. summary(s.reg2)
     
142. 
     
143. 
     
144. 
     
145. 
     
146. 
     
147. ## multilevel logits
     
148. 
     
149. library(foreign)
     
150. votedata <- read.dta("polls.dta")
     
151. 
     
152. pooled.logit1 <- glm(bush ~ black + female + edu + age, family=binomial(link="logit"), data=votedata)
     
153. summary(pooled.logit1)
     
154. 
     
155. 
     
156. ## Multilevel logit -- random intercept
     
157. multi.logit1 <- glmer(bush ~ black + female + edu + age + (1 | state), family=binomial(link="logit"), data=votedata)
     
158. summary(multi.logit1)
     
159. 
     
160. 
     
161. ## Multilevel logit -- random slopes
     
162. 
     
163. multi.logit2 <- glmer(bush ~ black + female + age + (1 + edu | state), family=binomial(link="logit"), data=votedata)
     
164. summary(multi.logit2)
     
165. 
     
166. 
     
167. ## Using Zelig
     
168. 
     
169. # logit example 1
     
170. 
     
171. z.logit1 <- zelig(bush ~ black + female + edu + age + tag(1 | state), data=votedata, model="logit.mixed")
     
172. summary(z.logit1)
     
173. 
     
174. x.high <- setx(z.logit1, age=4)
     
175. x.low <- setx(z.logit1, age=1)
     
176. 
     
177. s.logit1 <- sim(z.logit1, x=x.high, x1=x.low)
     
178. summary(s.logit1)
     
179. 
     
180. # logit example 2 -- note assumed intercept in random component
     
181. 
     
182. z.logit2 <- zelig(bush ~ black + female + age + tag(edu| state), data=votedata, model="logit.mixed")
     
183. summary(z.logit2)
     
184. 
     
185. ## something goes wrong from here with setx command
     
186. 
     
187. #x.high <- setx(z.logit2, age=4)
     
188. #x.low <- setx(z.logit2, age=1)
     
189. 
     
190. #s.logit2 <- sim(z.logit2, x=x.high, x1=x.low)
     
191. #summary(s.logit2)

复制代码

1. ## multiple imputation code for R ##
   
2. 
   
3. ## Example 1 ##
   
4. 
   
5. library(foreign)
   
6. CCSData <- read.dta("CCS2010v2.dta", convert.factors=FALSE)
   
7. 
   
8. install.packages("Rcpp") # we need to do this because otherwise Amelia might load a defunct version that doesn't work
   
9. 
   
10. install.packages("Amelia")
    
11. library(Amelia)
    
12. 
    
13. ## imputation ##
    
14. 
    
15. impute.out <- amelia(CCSData, m=5, noms=c("develop1","develop2"))
    
16. 
    
17. # built-in diagnostics #
    
18. 
    
19. summary(impute.out)
    
20. 
    
21. missmap(impute.out)
    
22. 
    
23. plot(impute.out)
    
24. 
    
25. overimpute(impute.out, var="age")
    
26. 
    
27. ## save imputed datasets ##
    
28. 
    
29. write.amelia(obj=impute.out, file.stem = "impdata", format = "dta")
    
30. 
    
31. impd1 <- subset(impute.out$imputations[[1]])
    
32. impd2 <- subset(impute.out$imputations[[2]])
    
33. impd3 <- subset(impute.out$imputations[[3]])
    
34. impd4 <- subset(impute.out$imputations[[4]])
    
35. impd5 <- subset(impute.out$imputations[[5]])
    
36. 
    
37. 
    
38. install.packages("Zelig")
    
39. library(Zelig)
    
40. install.packages("ZeligChoice")
    
41. library(ZeligChoice)
    
42. 
    
43. #regmod1.out <- zelig(as.factor(develop1) ~ incgroup + rent + southcoast + NEP, model="mlogit", data=CCSData)
    
44. #summary(regmod1.out)
    
45. 
    
46. #impmod1.out <- zelig(as.factor(develop1) ~ incgroup + rent + southcoast + NEP, model="mlogit", data=mi(impd1, impd2, impd3, impd4, impd5))
    
47. #summary(impmod1.out)  # not working -- we can do it by hand
    
48. 
    
49. install.packages("nnet")
    
50. library(nnet)
    
51. 
    
52. regmod1.out <- summary(multinom(develop1 ~ incgroup + rent + southcoast + NEP, Hess=TRUE, data=CCSData))
    
53. coeffs <- cbind(t(regmod1.out$coefficients[1,]), t(regmod1.out$coefficients[2,]))
    
54. ses <- sqrt(diag(solve(regmod1.out$Hessian)))
    
55. zs <- coeffs/ses
    
56. ps <- 2*(1 - pnorm(abs(zs)))
    
57. 
    
58. reg.final <- t(rbind(coeffs, ses, zs, ps))
    
59. colnames(reg.final) <- c("Coeff", "SE", "Z", "P")
    
60. 
    
61. impmod1.out <- summary(multinom(develop1 ~ incgroup + rent + southcoast + NEP, Hess=TRUE, data=impd1))
    
62. impmod2.out <- summary(multinom(develop1 ~ incgroup + rent + southcoast + NEP, Hess=TRUE, data=impd2))
    
63. impmod3.out <- summary(multinom(develop1 ~ incgroup + rent + southcoast + NEP, Hess=TRUE, data=impd3))
    
64. impmod4.out <- summary(multinom(develop1 ~ incgroup + rent + southcoast + NEP, Hess=TRUE, data=impd4))
    
65. impmod5.out <- summary(multinom(develop1 ~ incgroup + rent + southcoast + NEP, Hess=TRUE, data=impd5))
    
66. 
    
67. coeffs1 <- cbind(t(impmod1.out$coefficients[1,]), t(impmod1.out$coefficients[2,]))
    
68. se1 <- sqrt(diag(solve(impmod1.out$Hessian)))
    
69. coeffs2 <- cbind(t(impmod2.out$coefficients[1,]), t(impmod2.out$coefficients[2,]))
    
70. se2 <- sqrt(diag(solve(impmod2.out$Hessian)))
    
71. coeffs3 <- cbind(t(impmod3.out$coefficients[1,]), t(impmod3.out$coefficients[2,]))
    
72. se3 <- sqrt(diag(solve(impmod3.out$Hessian)))
    
73. coeffs4 <- cbind(t(impmod4.out$coefficients[1,]), t(impmod4.out$coefficients[2,]))
    
74. se4 <- sqrt(diag(solve(impmod4.out$Hessian)))
    
75. coeffs5 <- cbind(t(impmod5.out$coefficients[1,]), t(impmod5.out$coefficients[2,]))
    
76. se5 <- sqrt(diag(solve(impmod5.out$Hessian)))
    
77. 
    
78. mi.coeffs <- rbind(coeffs1, coeffs2, coeffs3, coeffs4, coeffs5)
    
79. mi.se <- rbind(se1, se2, se3, se4, se5)
    
80. 
    
81. mi.results <- mi.meld(mi.coeffs, mi.se)
    
82. mi.results
    
83. 
    
84. mi.z <- mi.results$q.mi/mi.results$se.mi
    
85. mi.p <- 2*(1 - pnorm(abs(mi.z)))
    
86. 
    
87. mi.final <- cbind(t(mi.results$q.mi), t(mi.results$se.mi), t(mi.z), t(mi.p))
    
88. colnames(mi.final) <- c("Coeff", "SE", "Z", "P")
    
89. 
    
90. # without imputation
    
91. reg.final
    
92. #with imputation
    
93. mi.final
    
94. 
    
95. ## Example 2 ##
    
96. 
    
97. MIData2 <- read.dta("development.dta", convert.factors=FALSE)
    
98. 
    
99. MIData2sub <- subset(MIData2, select=c(gxpdhlth, gini, g, glag, dictator, births, infmort, cath, moslem, femsec, fertil,country,year))
    
100. 
     
101. impute2.out <- amelia(MIData2sub, m=5, cs="country", ts="year", noms="dictator")
     
102. 
     
103. tscsPlot(impute2.out, var="infmort", cs="1")
     
104. 
     
105. impute3.out <- amelia(MIData2sub, m=5, cs="country", ts="year", noms="dictator", polytime=3)
     
106. 
     
107. tscsPlot(impute3.out, var="infmort", cs="1")
     
108. 
     
109. regmod2.out <- zelig(infmort ~ gxpdhlth + glag + dictator + femsec, model="ls", data=MIData2)
     
110. summary(regmod2.out)
     
111. 
     
112. impd1 <- subset(impute3.out$imputations[[1]])
     
113. impd2 <- subset(impute3.out$imputations[[2]])
     
114. impd3 <- subset(impute3.out$imputations[[3]])
     
115. impd4 <- subset(impute3.out$imputations[[4]])
     
116. impd5 <- subset(impute3.out$imputations[[5]])
     
117. 
     
118. impmod2.out <- zelig(infmort ~ gxpdhlth + glag + dictator + femsec, model="ls", data=mi(impd1, impd2, impd3, impd4, impd5))
     
119. summary(impmod2.out)
     
120. 
     
121. impmod2.out <- zelig(g ~ gxpdhlth + glag + dictator + femsec, model="ls", data=mi(impd1, impd2, impd3, impd4, impd5))
     
122. summary(impmod2.out)

复制代码

1. ## Selection models in R ##
   
2. 
   
3. install.packages("sampleSelection")
   
4. library(sampleSelection)
   
5. 
   
6. data("Mroz87")
   
7. Mroz87$kids <- (Mroz87$kids5 + Mroz87$kids618 > 0)
   
8. head(Mroz87)
   
9. 
   
10. # Female labor supply (lfp = labour force participation)
    
11. 
    
12. ## Outcome equations without correcting for selection
    
13. # I() means "as-is" -- do calculation in parentheses then use as variable
    
14. 
    
15. ## Comparison of linear regression and selection model
    
16. 
    
17. outcome1 <- lm(wage ~ exper, data = Mroz87)
    
18. summary(outcome1)
    
19. 
    
20. selection1 <- selection(selection = lfp ~ age + I(age^2) + faminc + kids + educ, outcome = wage ~ exper,
    
21. data = Mroz87, method = "2step")
    
22. summary(selection1)
    
23. 
    
24. plot(Mroz87$wage ~ Mroz87$exper)
    
25. curve(outcome1$coeff[1] + outcome1$coeff[2]*x, col="black", lwd="2", add=TRUE)
    
26. curve(selection1$coeff[1] + selection1$coeff[2]*x, col="orange", lwd="2", add=TRUE)
    
27. 
    
28. 
    
29. ## A more complete model comparison
    
30. 
    
31. outcome2 <- lm(wage ~ exper + I( exper^2 ) + educ + city, data = Mroz87)
    
32. summary(outcome1)
    
33. 
    
34. ## Correcting for selection
    
35. 
    
36. selection.twostep2 <- selection(selection = lfp ~ age + I(age^2) + faminc + kids + educ, outcome = wage ~ exper + I(exper^2) + educ + city,
    
37. data = Mroz87, method = "2step")
    
38. summary(selection.twostep2)
    
39. 
    
40. selection.mle <- selection(selection = lfp ~ age + I(age^2) + faminc + kids + educ, outcome = wage ~ exper + I(exper^2) + educ + city,
    
41. data = Mroz87, method = "mle")
    
42. summary(selection.mle)
    
43. 
    
44. 
    
45. ## Heckman model selection "by hand" ##
    
46. 
    
47. seleqn1 <- glm(lfp ~ age + I(age^2) + faminc + kids + educ, family=binomial(link="probit"), data=Mroz87)
    
48. summary(seleqn1)
    
49. 
    
50. ## Calculate inverse Mills ratio by hand ##
    
51. 
    
52. Mroz87$IMR <- dnorm(seleqn1$linear.predictors)/pnorm(seleqn1$linear.predictors)
    
53. 
    
54. ## Outcome equation correcting for selection ##
    
55. 
    
56. outeqn1 <- lm(wage ~ exper + I(exper^2) + educ + city + IMR, data=Mroz87, subset=(lfp==1))
    
57. summary(outeqn1)
    
58. 
    
59. ## compare to selection package -- coefficients right, se's wrong
    
60. summary(selection.twostep2)
    
61. 
    
62. 
    
63. ## interpretation
    
64. ## If our independent variables does not appear in the selection equation, we can interpret beta as in linear regression
    
65. ## If it does appear in the selection equation, we must calculate:
    
66. 
    
67. beta.educ.sel <- selection.twostep2$coefficients[6]
    
68. beta.educ.out <- selection.twostep2$coefficients[10]
    
69. beta.IMR <- selection.twostep2$coefficients[12]
    
70. delta <- selection.twostep2$imrDelta
    
71. 
    
72. marginal.effect <- beta.educ.out - beta.educ.sel * beta.IMR * delta
    
73. mr2 <- marginal.effect * Mroz87$educ
    
74. 
    
75. 
    
76. plot(Mroz87$wage ~ Mroz87$educ)
    
77. lines(mr2 ~ Mroz87$educ, type="l", col="green", lwd="2")
    
78. 
    
79. ## Selection with a binary outcome variable
    
80. ## Data from Kimball (2006)
    
81. 
    
82. library(foreign)
    
83. conflict.data <- read.dta("MissingLink_JPRfinal.dta",  convert.factors=FALSE)
    
84. conflict.data <- na.omit(conflict.data)
    
85. head(conflict.data)
    
86. 
    
87. ## A probit model for conflict
    
88. 
    
89. probit.model <- glm(conflict ~ relcap+contig+jtdem+jaut+pwrs+allform2, family=binomial(link=probit), data=conflict.data)
    
90. summary(probit.model)
    
91. 
    
92. install.packages("Zelig")
    
93. library(Zelig)
    
94. install.packages("ZeligChoice")
    
95. library(ZeligChoice)
    
96. 
    
97. 
    
98. selection.formula <- list(mu1 = conflict ~ relcap+contig+jtdem+jaut+pwrs,
    
99. mu2 = allform2 ~ relcap+logdist+contig+jtdem+jaut+sharerival)
    
100. 
     
101. selection.binary <- zelig(selection.formula, model = "bprobit", data = conflict.data)
     
102. summary(selection.binary)
     
103. 
     
104. x.contig <- setx(selection.binary, contig=1)
     
105. sim.binary1 <- sim(selection.binary, x = x.contig)
     
106. summary(sim.binary1)
     
107. plot(sim.binary1)
     
108. 
     
109. x.noncontig <- setx(selection.binary, contig=0)
     
110. sim.binary2 <- sim(selection.binary, x = x.contig, x1=x.noncontig)
     
111. summary(sim.binary2)
     
112. plot(sim.binary2)

复制代码

1. ## Models for BTSCS, Class 5
   
2. 
   
3. IRdata <- read.table("OR.txt", header=T, sep="\t")
   
4. 
   
5. IRdata <- as.data.frame(na.omit(IRdata))
   
6. 
   
7. View(IRdata)
   
8. 
   
9. ## Pooled logit ##
   
10. 
    
11. logitmod1 <- glm(dispute ~ dem + growth + allies + contig + capratio + trade, family=binomial(link="logit"), data=IRdata, x=TRUE)
    
12. summary(logitmod1)
    
13. 
    
14. ## Fixed effects logit ##
    
15. ## the pglm package can supposedly also do this, but it seems to crash ##
    
16. 
    
17. install.packages("glmmML")
    
18. library(glmmML)
    
19. 
    
20. logitfemod1 <- glmmboot(dispute ~ dem + growth + allies + contig + capratio + trade, family=binomial(link="logit"), data=IRdata, cluster=ordyid)
    
21. summary(logitfemod1)
    
22. 
    
23. ## list first 10 estimated fixed effects ##
    
24. 
    
25. logitfemod1$frail[1:10]
    
26. 
    
27. ## how many observations in each model?
    
28. 
    
29. length(logitmod1$fitted.values)
    
30. 
    
31. dyadrep <- as.vector(table(IRdata$ordyid))
    
32. validobs <- as.numeric(logitfemod1$frail > -Inf)
    
33. sum(dyadrep * validobs)
    
34. 
    
35. ## dynamics
    
36. ## Logit with time dummy variables ##
    
37. 
    
38. logitmod2 <- glm(dispute ~ dem + growth + allies + contig + capratio + trade + as.factor(py), family=binomial(link="logit"), data=IRdata)
    
39. summary(logitmod2)
    
40. 
    
41. ## Logit with time effects: cubic polynomial ##
    
42. 
    
43. IRdata$py2 <- (IRdata$py)^2
    
44. IRdata$py3 <- (IRdata$py)^3
    
45. 
    
46. logitmod3 <- glm(dispute ~ dem + growth + allies + contig + capratio + trade + py + py2 + py3, family=binomial(link="logit"), data=IRdata)
    
47. summary(logitmod3)
    
48. 
    
49. ## plot hazard rate over time
    
50. 
    
51. xb <- (colMeans(cbind(1, IRdata$dem, IRdata$growth, IRdata$allies, IRdata$contig, IRdata$capratio, IRdata$trade))) %*% logitmod3$coefficients[1:7]
    
52. time <- seq(1,50,1)
    
53. 
    
54. hazard <- plogis(xb + (time*logitmod3$coefficients[8]) + ((time^2)*logitmod3$coefficients[9]) + ((time^3)*logitmod3$coefficients[10]))
    
55. plot(hazard~time, type="l")

复制代码

1. ## survival models in R ##
   
2. 
   
3. library(foreign)
   
4. 
   
5. coalition.data <- read.dta("coalition.dta",  convert.factors=FALSE)
   
6. coalition.data$fractionalization <- coalition.data$fractionalization/1000
   
7. 
   
8. ## duration measured in months
   
9. 
   
10. olsmodel1 <- lm(duration ~ investiture + fractionalization + polarization + majority_government + crisis, data=coalition.data)
    
11. summary(olsmodel1)
    
12. 
    
13. install.packages("Zelig")
    
14. library(Zelig)
    
15. 
    
16. ## the exponential survival model
    
17. 
    
18. exp.survival <- zelig(Surv(duration, censor12) ~ investiture + fractionalization + polarization + majority_government + crisis,
    
19. model="exp", data=coalition.data)
    
20. summary(exp.survival)
    
21. 
    
22. # expected values and first differences
    
23. x.minority <- setx(exp.survival, majority_government = 0)
    
24. x.majority <- setx(exp.survival, majority_government = 1)
    
25. 
    
26. exp.survival.sim <- sim(exp.survival, x=x.minority, x1=x.majority)
    
27. summary(exp.survival.sim)
    
28. 
    
29. plot(exp.survival.sim)
    
30. 
    
31. 
    
32. ## The Weibull model -- currently not working in Zelig, so we can work around it
    
33. 
    
34. #weibull.survival <- zelig(Surv(duration, censor12) ~ identifiability + volatility + response + investiture + polarization + fractionalization + majority_government, model="weibull", data=coalition.data)
    
35. 
    
36. weibull.survival <- survreg(Surv(duration, censor12) ~ investiture + fractionalization + polarization + majority_government + crisis,
    
37. dist="weibull", data=coalition.data, x=TRUE)  # note x=TRUE
    
38. summary(weibull.survival)
    
39. 
    
40. # scale is <1 so hazard rate decreasing over time
    
41. 
    
42. # plot hazard rate over time
    
43. p <- 1/weibull.survival$scale
    
44. t <- seq(1,60,1)
    
45. lambda <- exp(colMeans(weibull.survival$x) %*% weibull.survival$coefficients)
    
46. weibull.hazard <- lambda * p * (lambda * t)^(1-p)
    
47. plot(weibull.hazard, type="l")
    
48. 
    
49. 
    
50. # expected values and first differences
    
51. 
    
52. x.minority <- append(colMeans(weibull.survival$x),1)  # The 1 is so we can pick up the scale parameter
    
53. x.minority["majority_government"] <- 0
    
54. x.majority <- x.minority
    
55. x.majority["majority_government"] <- 1
    
56. 
    
57. weibull.coeffs <- append(weibull.survival$coefficients, log(weibull.survival$scale))
    
58. 
    
59. betas <- mvrnorm(1000, weibull.coeffs, vcov(weibull.survival))
    
60. 
    
61. expect.min <- 1/exp(-x.minority %*% t(betas))
    
62. expect.maj <- 1/exp(-x.majority %*% t(betas))
    
63. 
    
64. weibull.fd <- expect.maj - expect.min
    
65. 
    
66. mean.min <- mean(expect.min)
    
67. mean.maj <- mean(expect.maj)
    
68. sd.min <- apply(expect.min,1,sd)
    
69. sd.maj <- apply(expect.maj,1,sd)
    
70. mean.fd <- mean(weibull.fd)
    
71. sd.fd <- apply(weibull.fd,1,sd)
    
72. 
    
73. fd.results <- rbind(cbind(mean.min, sd.min), cbind(mean.maj, sd.maj), cbind(mean.fd, sd.fd))
    
74. colnames(fd.results) <- c("Mean", "SD")
    
75. rownames(fd.results) <- c("Minority Govt", "Majority Govt", "FD")
    
76. print(fd.results)
    
77. 
    
78. 
    
79. ## The Cox proportional hazard model -- once again not working in Zelig
    
80. 
    
81. #coxph.survival <- exp.survival <- zelig(Surv(duration, censor12) ~ identifiability + volatility + response + investiture + polarization + fractionalization + majority_government, model="coxph", data=coalition.data)
    
82. 
    
83. #install.packages("survival") -- already installed with Zelig
    
84. #library(survival)
    
85. 
    
86. coxph.survival <- coxph(Surv(duration, censor12) ~ identifiability + volatility + response + investiture + polarization + fractionalization + majority_government, data=coalition.data, x=TRUE)  # note x=TRUE
    
87. summary(coxph.survival)  # note reversal of signs!
    
88. 
    
89. # estimated survival function
    
90. 
    
91. plot(survfit(coxph.survival), xlab="Months", ylab="Governments Surviving")
    
92. 
    
93. x.minority <- colMeans(coxph.survival$x)  
    
94. x.minority["majority_government"] <- 0
    
95. x.majority <- x.minority
    
96. x.majority["majority_government"] <- 1
    
97. hyp.data <- as.data.frame(rbind(x.minority,x.majority))
    
98. 
    
99. plot(survfit(coxph.survival, newdata=hyp.data), xlab="Months", ylab="Governments Surviving", conf.int=T, lty=c(1,2))
    
100. legend("topright", legend=c("minority govt","majority govt"), lty=c(1,2))
     
101. 
     
102. # interpreting in terms of expected time to failure is hard because we do not estimate a baseline
     
103. # instead look at changes in the hazard rate
     
104. 
     
105. hazard.pct.change <- ((exp(x.minority %*% t(betas)) - exp(x.majority %*% t(betas)))/exp(x.minority %*% t(betas)))*100
     
106. hazard.mean <- mean(hazard.pct.change)
     
107. hazard.sd <- apply(hazard.pct.change,1,sd)
     
108. cbind(hazard.mean, hazard.sd)  # minority gov't 47% more likely to fail

复制代码

1. ## 2SLS in R ##
   
2. 
   
3. library(foreign)
   
4. 
   
5. AllData <- read.dta("newdail2002v2.dta",  convert.factors=FALSE)
   
6. CampaignData <- na.omit(subset(AllData, select=c(votes1st,incumb,spend_regular,spend_regularXinc,spend_public,partyquota1997,electoraK,dublin,senator,councillor,wonseat)))
   
7. 
   
8. ## OLS ##
   
9. 
   
10. olsmodel <- lm(votes1st ~ incumb + spend_regular + spend_regularXinc, data=CampaignData)
    
11. summary(olsmodel)
    
12. 
    
13. olsmodel2 <- lm(votes1st ~ incumb + spend_regular + spend_regularXinc + spend_public, data=CampaignData)
    
14. summary(olsmodel2)
    
15. 
    
16. ## 2SLS by hand -- this replicates models 1 and 2, Table 3 in Benoit and Laver 2008
    
17. 
    
18. first.stage.1 <- lm(spend_regular ~ partyquota1997 + electoraK + dublin + senator + councillor, data=CampaignData)
    
19. summary(first.stage.1)
    
20. 
    
21. CampaignData$instrumented.spending <- first.stage.1$fitted.values
    
22. CampaignData$inst.spend.inc <- CampaignData$instrumented.spending * CampaignData$incumb
    
23. 
    
24. second.stage.1 <- lm(votes1st ~ incumb + instrumented.spending + inst.spend.inc, data=CampaignData)
    
25. summary(second.stage.1)
    
26. 
    
27. 
    
28. first.stage.2 <- lm(spend_public ~ partyquota1997 + electoraK + dublin + senator + councillor, data=CampaignData)
    
29. summary(first.stage.2)
    
30. 
    
31. CampaignData$instrumented.pspending <- first.stage.2$fitted.values
    
32. 
    
33. second.stage.2 <- lm(votes1st ~ incumb + instrumented.spending + inst.spend.inc + instrumented.pspending, data=CampaignData)
    
34. summary(second.stage.2)
    
35. 
    
36. 
    
37. ## Issue to consider -- is interaction also endogenous?
    
38. ## Example of difference in standard errors
    
39. 
    
40. second.stage.1v2 <- lm(votes1st ~ instrumented.spending, data=CampaignData)
    
41. summary(second.stage.1v2)
    
42. 
    
43. ## 2SLS ##
    
44. ## compare standard errors
    
45. 
    
46. install.packages("AER")
    
47. library(AER)
    
48. 
    
49. tslsmodel <- ivreg(votes1st  ~ spend_regular | partyquota1997 + electoraK + dublin + senator + councillor, data=CampaignData)
    
50. summary(tslsmodel)
    
51. 
    
52. 
    
53. 
    
54. 
    
55. 
    
56. ## Complete models with 2SLS
    
57. ## note incumbency is also an instrument
    
58. 
    
59. tslsmodel.1 <- ivreg(votes1st  ~ incumb + spend_regular + spend_regularXinc | incumb + partyquota1997 + electoraK + dublin + senator + councillor, data=CampaignData)
    
60. summary(tslsmodel.1)
    
61. 
    
62. tslsmodel.2 <- ivreg(votes1st  ~ incumb + spend_regular + spend_regularXinc + spend_public | incumb + partyquota1997 + electoraK + dublin + senator + councillor, data=CampaignData)
    
63. summary(tslsmodel.2)
    
64. 
    
65. 
    
66. ## An alternate package
    
67. 
    
68. install.packages("sem")
    
69. library(sem)
    
70. 
    
71. tslsmodel.3 <- tsls(votes1st  ~ incumb + spend_regular + spend_regularXinc, ~ incumb + partyquota1997 + electoraK + dublin + senator + councillor, data=CampaignData)
    
72. summary(tslsmodel.3)
    
73. 
    
74. tslsmodel.4 <- tsls(votes1st  ~ incumb + spend_regular + spend_regularXinc + spend_public , ~ incumb + partyquota1997 + electoraK + dublin + senator + councillor, data=CampaignData)
    
75. summary(tslsmodel.4)
    
76. 
    
77. 
    
78. ## 2 stage estimator with a binary dv
    
79. 
    
80. probit.model <- glm(wonseat ~ incumb + spend_regular + spend_regularXinc + spend_public, family=binomial(link=probit), data=CampaignData)
    
81. summary(probit.model)
    
82. 
    
83. 
    
84. ## (sort of) replication of table 5
    
85. 
    
86. tspmodel.1 <- glm(wonseat ~ incumb + instrumented.spending + inst.spend.inc + instrumented.pspending, family=binomial(link=probit), data=CampaignData)
    
87. summary(tspmodel.1)
    
88. 
    
89. ## 2SCML
    
90. 
    
91. reg.inst1 <- lm(spend_regular ~ incumb + partyquota1997 + electoraK + dublin + senator + councillor, data=CampaignData)
    
92. reg.inst2 <- lm(spend_regularXinc ~ incumb + partyquota1997 + electoraK + dublin + senator + councillor, data=CampaignData)
    
93. reg.inst3 <- lm(spend_public ~ incumb + partyquota1997 + electoraK + dublin + senator + councillor, data=CampaignData)
    
94. 
    
95. CampaignData$inst1 <- reg.inst1$residuals
    
96. CampaignData$inst2 <- reg.inst2$residuals
    
97. CampaignData$inst3 <- reg.inst3$residuals
    
98. 
    
99. tscml.model <- glm(wonseat ~ incumb + spend_regular + spend_regularXinc + spend_public + inst1 + inst2 + inst3, family=binomial(link=probit), data=CampaignData)
    
100. summary(tscml.model)
     
101. 
     
102. ## test for endogeneity
     
103. 
     
104. install.packages("lmtest")
     
105. library(lmtest)
     
106. 
     
107. lrtest(probit.model, tscml.model)

复制代码